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- - - ABSTRACT 
0\ ' 

■ We obtain renormalized sets of right and left eigenvectors of the flux vector Jacobians of the rela- 

I tivistic MHD equations, which are regular and span a complete basis in any physical state including 

, degenerate ones. The renormalization procedure relies on the characterization of the degeneracy 

, ' types in terms of the normal and tangential components of the magnetic field to the wavefront in 

^ ' the fluid rest frame. Proper expressions of the renormalized eigenvectors in conserved variables are 

obtained through the corresponding matrix transformations. Our work completes previous analysis 
that present different sets of right eigenvectors for non-degenerate and degenerate states, and can be 
' seen as a relativistic generalization of earlier work performed in classical MHD. Based on the full wave 

fN| , decomposition (FWD) provided by the the renormalized set of eigenvectors in conserved variables, we 

have also developed a linearized (Roe-type) Riemann solver. Extensive testing against one- and two- 
dimensional standard numerical problems allows us to conclude that our solver is very robust. When 
compared with a family of simpler solvers that avoid the knowledge of the full characteristic structure 
of the equations in the computation of the numerical fluxes, our solver turns out to be less diffusive 
than HLL and HLLC, and comparable in accuracy to the HLLD solver. The amount of operations 
needed by the FWD solver makes it less efficient computationally than those of the HLL family in 
one-dimensional problems. However its relative efficiency increases in multidimensional simulations. 
O ' Subject headings: MHD - relativity - methods: numerical 

+i ■ 

C/3 
. 

_1j 1. INTRODUCTION 

Astrophysical plasma is probably the most exciting, complex, physically entangled state of matter in Nature. Even 
, under the restrictive assumption that such plasma can be modeled as a fluid, rather than as a statistical ensemble of 
. interacting particles, the diversity of physical phenomena it can host deserves a careful study. This is particularly true 
if the magnetic flelds providing the coupling for the plasma constituents become dynamically important. Additional 
_ challenges for the physical modeling of astrophysical plasma pop up if its constituents are relativistic, if the plasma 
. moves either at bulk speeds very close to the speed of light, or if it is embedded in ultra-intense gravitational fields. 
• ' This is so because relativistic flows in association with intense gravitational and magnetic fields are commonly linked 
up to extremely energet ic phenomena in the Universe, viz. pulsar winds |Rces & Gunn 1974), anomalous X-ray pulsars 
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jyan Para diis et "all 1995), soft gamma-ray repeaters (SGR; 'Kouveliotou et al.i ll998l) . gamma-ray bursts (GRBs; E 

iranl 

12004 ) ■ relativistic jets in active galactic nuclei (AGNs; Bcgelman ct al. 1984j), etc. 

Neutron stars are among the astrophysical objects, which may harbor ultra-intense magn etic fields. Dipole magnetic 
fields in exces s of 10^^ G have been inferred from the spin-down rates of radio pulsars (jTavlor fc Stinebrind 119861 : 
lKulkarni|[l99l ). As the total magnetic energy of the star is probably much greater than the exterior dipole component, 
it is expected that the neutron star interior will be pierced by even larger field values 10^"* G). Also huge magnetic 
fields {'-^ 10^'* G) are attributed to anomalous X-ray pulsars (jKouveliotou et aI..1998., l. Extreme magnetars as SGRs host 
dipole magnetic fields as large as 10^* — 10^^ G, and interior fields in the range (5 - 10) X 10^5 G, ([Thompson fc Duncanl 
119961 : IWoods fc ThompsonI 120061 ). The mechanism to build up such superlative fields is still disputed, but they shall 
be produced in the course of the collapse of massive stellar cores, not only as a result of the compression of the stellar 
progenitor field, but also, most likely , because of the e fficient field amplification by means of the magneto-rotational 
instability fe.g.. [ Akivama et al.l [20031 : lObergaulinger et al. 2009b) or the dynamo action of in the proto neutron star 
([Thompson fc Duncanlll993D . 

Low-mass X-ray binaries containing a black-hole candidate display high-frequency quasi-periodic oscillations (QPOs) 
in their spectra. The QPO activity has been theoret ically attributed to t he normal modes of oscilla t ion of thick 
accretion disks (tori) girding stellar-mas black holes ([Rezzolla et al.l I2003[) . Recently, IMontero et al.l ([2007D have 
shown that such oscillations also happen if the tori (assumed to obey an ideal gas equation of state) are pierced by 
toroidal magnetic fields of non - negligible strength; the f undamental mode and its overtones being very similar to those 
of unmagnetized tori (but see lBarkov fc BaushevI 120091 . for the damping properties that a realistic equation of estate 
and more elaborated microphysics may have on the normal modes of oscillation). Q POs have also been observed in 
the X-ray tail of giant fiares of SGRs ([Israel et al.l 120051 : iStrohmaver fc Watts! I200"5f l. These flares are attributed to 
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crust-quakes in magnetars (|Duncan fc Thompsoiilfl996[ ) . Recent general relativistic m agnetohydrodynamic (GRMHD) 
studies show that the QPOs can be created by the oscillation mode s of the crust (jSamuelsson fc AnderssonI l2007f) 
coupled to Alfven oscillations of the interior (jCerda-Duran et"aLll2009f) . 

So me models o f GRB progenitors link them with the formation of rapidly rotating, highly magnetized neutron 
stars ()Usovlll99^ . Rotational and/or magnetic energy (a few 10^^ erg), decimated in scales of seconds, in a pulsar-like 
mechanism is released in the form of an initially opaque fireball. The 7— radiation appears via plasm a instabilities when 
the MHD approximation of the pulsar wind breaks down at scales of ~ 10^^ cm ()Thompsonlll994[ ). Alternatively, the 
rotational energy of a newly-born stellar-mass black hole c an be tap into the plasma above the poles of the black hole by 
gigantic magnetic fields threading its ergosph ere ([Bl andfor d fc Znaieklll977f). This mechani s m has been proven to work 
on a number of GRMHD simulations (e.g., iMcKinnevI 120061 : iBarkov fc Komissarovl 120081 . 120091 : iTchekhovskov et all 
|2008| ) provided that both the s pecific angular momeiitum of the black hole, and the magnetic field strength and 
topology are adequate (see also lAlov fc Obergaulingerl 120071). For the subclass of short-duration GRBs a merger of 
compact objects in a close binary is the progenitor of choice ( Paczvns"killl986l : lGoodmaiilll986l : [Mochkovitch et a"nil993l . 
, Aloy, Janka fc Miiller 2005). Very recently, the role played by the magnetic field in the merger of neutron stars has 
started to be invest igated b y means of global GRM HD simulations ()Price fc Rosswogjl2006l : lRosswosJl2007t iLiu et all 
l2008t [Anderson et~a l. 2008; Giacomazzo et al.) 120091 ). These simulations suggest that intense magnetic fields, in some 
cases exceeding by far 10^^ G, may develop because of the Kelvin-Helmholtz instabilities rising form the contact surface 
between the two neutron stars. However, the resolution of the afore mentioned global ap proaches is too coarse to extract 
definitive conclusions, as demonstrated by local MHD simulations (jObergaulinger et a l. 2009a). 

Magnetic fields might also be non-negligible, even at large distances from the formation site, in relativistic outflows 
associated to pulsar winds, AGNs, microquasars and GRBs. Even if the MHD app roximation may break down 
when modeling pulsar magnetospheres, it is a good starting point for their study (|Komissarovl 120061 ). This has 
motivated a number of works where relativistic magneto-hydrodynamic (RMHD) simulations have been employed as 
tools to increase our understanding of the flow structure, observed emission, polarization and spectral properties (e.g., 
iBucciantini et aLll2006l : lKomissarovll2006l ). The plasma flowing along the channels drilled by relativistic jets in AGNs, 
microquasars and GRBs should be magnetized to some extent, since we observe non-thermal spectra that can be 
partly produced by the synchrotron radiation emitted by leptons gyrating in a magnetic field (see, e.g.. lFerrarilll998l 
for a review). The high degree of polarization of many AGNs at parsec scale (e.g., [Gabuzda et al. 200(|) advocates in 
favor of magnetic fields not being (only) randomly orient ed. At larger scales linear polariz ation maps reveal a certain 
ordered topology, which vari es from source to sourc e (e.g., IWalker et al.lll987l : rKigure et al.ir2004 etc.). For the galactic 
microquasar GRS1915-hl05,[RodrI guez fc Mirabel (| 1999( 1 find magnetic fields ~ 0.05 G, to account for the observed 
spectrum. Likewise, the ejecta associated to GRBs might be moderately or highly magnetized. The paucity of optical 
flashes in the afterglow of most GRBs can be readil y explained if the magnetization of the flow (measured by the 
Foynting to inertial flux ratio) is close to unity (e.g., Giannios et al.|[2b08l : iMimica "eran i200^. 

The necessity to model the aforementioned astrophysical scenarios in the framework of (general) RMHD, together 
with the fast increase in computing power, is pushing towards the development of more and more accurate and efficient 
numerical algorithms. In the last years, considerable progress has been achieved in numerical special RMHD, mainly 
trough the extension of high-resolution shock- capturing (HRSC) me thods, after the success of this kind of techniques 
in special relativistic hydrodynamics (see, e.g., lMartf fc Miilleill2003l ). HRSC methods are written in conservation form 
and the time evolution of zone averaged state variables is governed by some functions, the numerical fluxes, evaluated at 
zone interfaces. In the so called Godunov-type methods, an important subsample of HRSC methods, numerical fluxes 
are evaluated through the exact or approximate solution of the Riemann problem, the decay of an initia l discontinu ity 
between constant states. Despite the fact that such an analytical solution in RMHD is known (Romero et al.ll2605l for 
some particular orientation of the magnetic fleld and the fluid velocity fleld; [Giacomazzo fc Rezzolla 2006), approximate 
algorithms are usually preferred (because of their larger numerical efficiency). Komissarov (1999), Balsara (2001), and 
iKoldoba et al. (200^ developed independent algorithms {Roe-type algorithms; sec Ma rti fc Miillerl 120031 ) based on 
linearized Riem ann solvers relyin g on th e characteristic structure of the RMHD equations. 

More recently, iDel Zanna et al.l ()2003[ ) have developed a third order HRSC scheme for RMHD which does not require 
the knowledg e of the characteristic structure of the equations to obtain the numerical fluxes. The algorithm follows the 
iHarten et al.l (1983; HLL) approach that approximates the full structure of the Riemann solution (with seven waves 
separating the two initial states plus six intermediate states ) by just two li miting waves with a single intermediate state. 
The algorithms developed by Leismann et al. (2005) and v an der Hoist et al. (2008) are also b ased in the simple HLL 
appro ach. In an effort to improve the accuracy of the HLL strategy. iMignone fc Bodol (|2009) and lHonkkila fc Janhuneii 
( 2007 ) have extended to RM HD the so called HLLC scheme (initially devised for the Euler equations by iToro et al. 
(| 19941 ). and extended by[Gurski (2 00l andlO (|2005f l to classical MHD). In the HLLC scheme, the middle (entropy) 
wave is also captured. Finally, IMignone et al.l (2009; MUB09) have developed a flve-wave HLL solver for RMHD to 
capture Alfven discontinuities. 

Following the trail of special relativistic MHD, numerical GRMHD is blossoming. Komissarov' (200 5i) has e xtended 
his code based in the linearized Riemann s o lver to GRMHD us i ng the approach described in Pons et ^ ahl (Il998,). 
The codes developed bv lGammie et al.l ()2003[ ), lMizuno et all (|2006l) . ITchekhovskov et al.1 (I2007D. and Del Z anna et al] 
(pool are based in the HLL approach. The interested reader is addressed to the review of lFontI (2008. ) for an up-to- 
date list and description of codes in both general relativistic hydrodynamics and magnet ohydr odynamics, including 
approaches not mentioned in this introduction, as symmetric TVD schemes (jKoide et al.lll999| ) or artificial viscosity 
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methods iDeVilliers fc Ha^devl (|200l . 

The purpose of the present paper is twofold. On one hand, to present a renormalized set of right and left eigenvectors 
of the flux vector Jacobians of the RMHD equations which are regular and span a complete basis in any physical state, 
including degenerate states. On the other hand, to evaluate numerically the performance of a RMHD Riemann solver 
based on the aforementioned spectral decomposition. Both the theoretical analysis and the numerical applications 
presented in this paper are based on the work developed by IAnt6n| (i20081. Moreov er, the numerical procedure 
d escribed in this pa per was already used in GRMHD calculations bv I Anton etaH (I2OO6I) using the method described 
in lPons et al.l ()1998[ ) to extend the RMHD solver to GRMHD. Our nu merical approach deviates in several aspects from 
previ ous works based on linearized Riemann solver approaches (i.e., iKomissarovl 1 1999t iBalsaral [20011 : iKoldoba et al.l 
I2OO2I ). First, numerical fluxes are computed from the spectral decomposition in conserved variables. Seco nd, we 
presen t explicit expressions also for the left eigenvectors. Third, and more important, we have extended Brio fc Wul 
([1988) classical MHD strategy to relativistic flows, giving sets of right and left eigenvectors which are well defined 
through degenerate states. 

The organization of the paper is as follows. In § 2 the equations of RMHD are presented and the suitable definitions 
of variables given. In § 3 the mathematical structure of the equations is reviewed. Special attention in this section 
is paid to the characterization of degeneracies and the procedure of renormalization of the right and left eigenvectors 
in covariant variables. Sect. 4 includes a brief note concerning the non-convex character of relativistic MHD. Sects. 5 
and 6 are devoted to obtain, respectively, the renormalized right and left eigenvectors in conserved variables. Sect. 7 
summarizes the properties of the numerical code based in the full wave decomposition Riemann solver built from the 
renormalized eigenvectors presented in the preceding sections. A thorough study of the performance of the code, based 
on a battery of one and two-dimensional tests and the comparison with other numerical strategies, is presented in § 8. 
Sect. 9 presents the conclusions of our work. 

2. THE EQUATIONS OF IDEAL RELATIVISTIC MAGNETOHYDRODYNAMICS 

Let J^, r^'' and *F^^^ {n, v = 0,1, 2, 3) be the components of the rest-mass current density, the energy-momentum 
tensor and the Maxwell dual tensor of an ideal magneto-fluid, respectively 

= pu" (1) 
Tf"" = ph*ut'u'' + gf^p* - b^b" (2) 

*Ff' =u''f' -u^b^", (3) 

where p is the proper rest-mass density, h* = 1 + e + p/ p + b"^ / p is the specific enthalpy including the contribution 
from the magnetic field (6^ stands for b^bfj,), e is the specific internal energy, p the thermal pressure, p* ^ p + 6^/2 the 
total pressure, and g^'^ the metric of the space-time where the fluid evolves. Throughout the paper we use units in 
which the speed of light is c = 1 and the (47r)^/^ factor is absorbed in the definition of the magnetic field. 

The four- vectors representing the fluid velocity, u^, and the magnetic field measured in the fluid rest frame, b'^, satisfy 
the conditions u''u^ = — 1 and u^^b^ = 0, and there is an equation of state relating the thermodynamics variables, p, p 
and e, p = p{p, e). All the discussion will be valid for a general equation of state but results will be shown for an ideal 
gas, for which p = (7 — where 7 is the adiabatic exponent. 

The equations of ideal RMHD correspond to the conservation of rest-mass and energy-momentum, and the Maxwell 
equations. In a flat space-time and Cartesian coordinates, these equations read: 

= (4) 

T^-;^ = (5) 

*F>'\^ = 0, (6) 

where subscript denotes partial derivative with respect to the corresponding coordinate, {t, x, y, z), and the standard 
Einstein sum convention is assumed. Greek indices will run from to 3 (or from t to z) while Roman run from 1 to 3 
(or from a; to 2;). 

The above system can be written as a system of conservation laws as follows 
where the state vector, U, and the fluxes, F* (j = 1, 2, 3 or i = x, y, z), are the following column vectors, 
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(8) 



TV' +p*v' -b"B'/W ■ ^ ' 

\ v'B'' - v'^B' / 

In the preceding equations, D, and r stand, respectively for the rest-mass density, the momentum density of the 
magnetized fluid in the j-direction and its total energy density as measured in the laboratory (i.e., Eulerian) frame, 

D = pW, (10) 
^ ph*W^v^ ^b%^, (11) 

T ^ ph*W^ -p* - {b^f - D. (12) 

Quantities f ' stand for the components of the fluid velocity trivector as measured in the laboratory frame. They are 
related with the components of the fluid four-velocity according to the following expression = W{l,v^ jV"^ ,v^), 
where W is the flow Lorentz factor, = 1/(1 — v'vi). 

The following fundamental relations hold between the components of the magnetic fleld four-vector in the comoving 
frame and the three vector components -B' measured in the laboratory frame, 



b° = WB-v, (13) 

= (14) 

V and B being, respectively, the tri- vectors (u^,ti^,i;^) and (B^ , B^ , B^). 
Finally, the square of the modulus of the magnetic field can be written as 

b^ = ^ + iB..r. (15) 

The preceding system must be complemented with the time component of equation ([S]), that becomes the usual 
divergence constraint 

dB' , , 

which should be fulfilled at all times. 

Fluxes F* (i — x, y, z) are functions of the conserved variables, U, although for the RMHD this dependence, in 
general, can not be expressed explicitly. It is therefore necessary to introduce another set of variables, the so-called 
primitive variables^ derived from the conserved ones, in terms of which the fluxes can be computed explicitly. The set 
of primitive variables used along this work is written as the column vector 

V={p,p,v^,vy,v^B'^,By,BY, (17) 

where the superscript stands for the transposition. The transformation between conserved and primitive variables is 
done from the following expressions (see Komissarov 1999, Leismann et al. 2005), obtained, after some algebra, from 
equations (fTTj) . (fT2)) . (fft jl and (fT4| . 

S' = {Z + B^r^^ {2Z + B2)^5_|)!, (18) 

where S is the tri-vector (5^, S"", S^) and Z = phW^ . 

Equations (|18p and (fT5|) . together with the definitions of D (equation [TU)) and Z, form a system for the unknowns 
p, p and assuming the function h = h(p,p) is provided. In our calculations, since we restrict ourselves to an ideal 
gas equation of state, h = 1 + 7p/p(7 — 1). 



3. CHARACTERISTIC STRUCTURE OF THE RMHD EQUATIONS 

The hyperbolicity of the equations of RMHD including the derivati on of wavespeeds and t he corresponding eigen- 
vector s, and the analysis of various degeneracies has been studied bv lAnile fc Pennisl (jl987| ) and reviewed bv lAnild 
(|1989f ). These authors did their analysis in a covariant framework, using a set of variables, the so-called covariant 
variables, in which, the vector of unknows. 



(20) 



is augmented to 10 variables, where s is the specific entro py. The discussion has been more recently remasterized and 
cast in a form more suitable for numerical applications bv lKomissarovl (jl999D (see also lBalsara|[200ll ). We shall review 

it here. 

In terms of variables U (see lAniiill989f) . the system of RMHD equations can be written as a quasi-linear system of 
the form 



y^^U;^ = 0, 

where stands for the covariant derivative, and the 10 x 10 matrices A'^ are given by 



(21) 



U^/c2 0^ 

0^ J 



phS'i 



0^ 



(22) 



where Cs stands for the speed of sound 



(23) 



e being the mass-energy density of the fluid e ~ p{l + e) . In equation p2p the following definitions are introduced: 



£ = ph + b^, 

= {phgf^^ + {ph - bycl)u^^u")/{ph), 
ft^°' = {u''bf'/cl~ut'b'^)/{ph), 



(24) 
(25) 
(26) 
(27) 



as well as the notation 



0^=0, 0"^ = (0,0,0,0)"^, 0^ = (0,0,0,0). 



(28) 



It is important to remark that the 10 covariant variables we have used to write the system of equations are not 
independent, since they are related by the constraints 



and 



u'^Ua — —1, 
= 0, 

ac«(u"6°-u*'6") =0, 



(29) 
(30) 
(31) 



which reduces to the usual divergence constraint. 



3.1. Wavespeeds 

If (j){x'^) — defines a characteristic hypersurface of the above system (HI]), the characteristic matrix, given by A°'4>a 
can be written as 



^SaS^^mi^ 0^^ 
_ I B6i; aSi^ 0^ 
^ \ ph(j), 0, a/cl 

0^ 0^ a 



(32) 
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where 0^ = 9^0, a = B = b°'(t>a, G = = l^"^(f>a = (f)'' + (ph-b^ /c^^)au^' /{ph) + Bbf" / {ph), = /^"</>„ = 

{ab'^/c^ — Bu^)/{ph), and to(^ = (0^ + 2au'^)bi, — S(5^. Since = is a characteristic surface, the determinant of 

the matrix (|32p must vanish, i.e. 

det(^"(/.„) = £ a^A'^Mi = , (33) 

where 



A = £a^-B^ 



1 



Afi^phl — -1] a 



ph+—] a'G + B'G. 



(34) 
(35) 



The above equations, vahd for a general space-time, will be applied to obtain the wavespeeds in a flat space-time in 
Cartesian coordinates. To this end, let us consider a planar wave propagating in a given direction, that we choose as 
the a;-axis, with speed A. The normal to the characteristic hypersurface describing this wave is given by the four-vector 



0^ = (-A, 1,0,0), (36) 

and by substituting equation p6p in equation (|33p we obtain the so called characteristic polynomial, whose zeroes give 
the characteristic speed of the waves propagating in the x-direction. Three different kinds of waves can be obtained 
according to which factor in equation p3p becomes zero. For entropic waves a = 0, for Alfven waves ^ = 0, and for 
magnetosonic waves Mi — 0. 

The characteristic speed A of the entropic waves propagating in the x-direction, given by the solution of the equation 
a = 0, is the following 



A = t;-(=Ae). (37) 
For Alfven waves, given by ^ = 0, there are two solutions corresponding, in general, to different speeds of the waves, 

A = ^^^^(=A„±). (38) 

In the case of magnetosonic waves there are four solutions, two of them corresponding to the slow magnetosonic waves 
and the other two to the fast magnetosonic waves. It is however not possible, in general, to obtain simple expressions 
for their speeds since they are given by the solutions of the quartic equation Mi — with a, B and G in equation (|35p 
explicitly written in terms of A as 



a = iy(-A + w^), 
i3 = 6^ -6"A, 
G = l- A^ 



(39) 
(40) 
(41) 



It is worth noticing that just as in the classical case, the seven eigenvalues corresponding to the entropic (1), Alfven 
(2), slow magnetosonic (2) and fast magnetosonic waves (2) can be ordered as follows 

A7 < A- < a; < Ae < A+ < A+ < A+, (42) 

where the subscripts e, a, s and / stand for entropic, Alfven, slow magnetosonic and fast magnetosonic respectively, 
and the superscript — or -I- refer to the lower or higher value of each pair. This fact allows us to group the Alfven and 
magnetosonic eigenvalues in two classes separated by the entropic eigenvalue. It should be also emphasized that the 
superscript does not correspond to the sign on the left hand side of equation ((38)) . We will use subscripts to refer to 
this sign whenever needed. 

Let us remark that in the previous discussion about the roots of the characteristic polynomial we have omitted 
the fact that the entropy waves as well as the Alfven waves appear as double roots. These superfluous eigenvalues 
are associ ated with unph ysical waves and are the result of working with the unconstrained system of equations. We 
note that Ivan PuttenI (jT991) derived a different augmented system of RMHD equations in constrained-free form with 
different unphysical waves. Any attempt to develop a numerical procedure to solve the RMHD equations based on 
the wave structure of the RM HD equations must remo ve these unphy s ical w aves (i.e. the corresponding eigenvectors) 
from the wave decomposition. IKomissarovl (|1999f ) and IKoldoba et al.l ([20021 ) eliminate the unphysical ei genvectors b y 
demanding the waves to preserve the values of the invariants u'^u^ = — 1 and = as suggested bv lAnild (|1989f l. 
Correspondingly, Balsara, (2001. ) selects the physical eigenvectors by comparing with the equivalent expressions in the 
nonrelativistic limit. 
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In Appendix [X] we discuss our results obtained in the representation of the characteristic speeds for fluids under 
different thermodynamical conditions and states of motion. These diagrams show the normal speed of planar wavefronts 
propagating in different directions (phase speed diagrams). These diagrams are the genera lization to relativistic, 
arbitr arily moving flows of the o riginal diagrams introduced by Friedrichs (see (jJeffrev fc Tan iuti 1964)). In a recent 
paper. [Keppens &: MeliaiiH ([2008) have revisited the theory for linear RMHD wave propagation showing the equivalence 
with the characteristic speed approach, and the connection between the phase and group speed diagrams by means of 
a Huygens construction, in arbitrary reference frames. 

3.2. Degeneracies 

As in the case of classical magnetohydrodynamics, the equations of RMHD present degeneracies in the sense that 
two or more eigenvalues of the system of equations coincide and, therefore, the strict hyperb olicity of the system breaks 
down. The conditions for that to happen have been analyzed in the context of RMHD by iKomissciro vl ([l999l) and we 
summarize them here. Komissarov's analysis is done for the characteristic wave speeds in the fluid frame and leads 
to the conclusion that the conditions of degeneracy in RMHD are the same as in classical MHD. Degeneracies are 
encountered for waves propagating perpendicular to the magnetic field direction (Type I) and for waves propagating 
along the magnetic field direction (Type II). Finally, a particular subcase of Type II degeneracy appears when the 

sound speed is equal to Ca = \fWl£. We will refer to this special subcase as Type II'. 

For Type I degeneracy, the two Alfven waves, the entropic wave and the two slow magnetosonic waves propagate 
at the same speed (A~ = Aj = Ae = A^ = A^). For Type II degeneracy, an Alfven wave and a magnetosonic wave 
(slow or fast) of the same class propagate at the same speed (Aj = A^ or A^ — or A+ = A+ or A+ = A^ ). In the 
special Type II' subcase, an Alfven wave and both the slow and fast magnetosonic waves of the same class propagate 
at the same speed {XJ = A^ = Aj or A+ = A+ = A^). We have to emphasize that in classical MHD, if Type II 
degeneracy appears, both Alfven waves exhibit such a degeneracy. In RMHD however, due to aberration, the condition 
for degeneracy can be fulfilled for one Alfven wave but not for the other. Only if the tangential component of the 
fluid velocity vanishes, we recover the classical behaviour. In Appendix [X] we present a series diagrams displaying the 
characteristic wave speeds as a function of direction for fluid with different thermodynamical conditions and states of 
motion, payin g atten tion to the manifestation of the various types of degeneracy in the laboratory frame. 

iKomissarovl ()1999D also gives a covariant characterization of the different types of degeneracy. Type I degeneracy is 
characterized by B, equation (j40p . being zero for the case of the Alfven wave, whereas Type II degeneracy arises when 
the quantity B'^/ (G+a^) for the Alfven wave equals to b^. We have made a step forward concerning the characterization 
of both types of degeneracy by introducing the components of the magnetic field parallel and perpendicular to the 
Alfven wave in the comoving frame. To this end, let v"' b e the unitary vector, normal to the front wave as described 
by the comoving observer with the ffuid fsee IAnil"ill989D . To obtain this vector, first we project the vector 0" onto 
the orthogonal space to the fluid four velocity and, second, we normalize this projection. It is easy to see that this 
vector is given by 

^ (43) 



which is well defined even for the degenerate cases. 

Now, we can characterize both types of degeneracies in terms of the magnetic field components normal and tangential 
to the Alfven wave front in the comoving frame, and bf respectively. The magnetic field 6" can be written as 

b'^=b"^+b?^ibp:^0Jiy: + b?, (44) 



where denotes the normal to the Alfven front wave. Using this decomposition and the definition of Va, equation (|43p 
particularized to the Alfven wave, it is easy to check that 

= « = (6.^,")' = ^^. (45) 

Now, Type I degeneracy reduces to the case when 6^ — (0,0,0,0), whereas Type II degeneracy is characterized by 
6^ = i.e.. Type I degeneracy occurs whenever the component of b" along the direction of propagation of the Alfven 
wave (normal to the Alfven wave front), 5^ is equal to zero, whereas Type II degeneracy occurs when the component 
of the magnetic field tangential to the Alfven wave front, bf, is equal to zero. Finally, let us note that the condition 
6" = (0, 0, 0, 0) leads, for a wave propagation along the x direction in the laboratory frame, to = 0. 

3.3. Renormalized right eigenvectors 

Let us start with the right eigenvectors for the eigenvalue problem {A'^(j)fj,)r = 0, as derived by Anile (1989), 
Entropy eigenvector. 

re = (0",0",0,1)T. (46) 

Alfven eigenvectors: 
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r,,± = (a£"^^,7.^0■'&^ Be%y<j,'>b', 0, 0)^, (47) 

where the dependence in the correspondmg wavespeeds (eigenvahies), Xa.±, is hidden in the quantities a, B and (jf , 
and the Levi-Civita ahernating tensor. 

Magnetosonic eigenvectors: 

r™,± = {ad'^,Bd'^ + aAr, c?A, 0)"^, (48) 

(m = s, /) where 



where, again, the dependence in the corresponding characteristic wavespeeds, As,± and A/^±, is hidden in the quantities 
a, B, G, A and Eigenvalues Xm.± are defined as belonging to the same class as Aa.±, i.e.. 



A, 



m,± 



if ^a,± = AJ. 



As it is well known, Alfven and magnetosonic eigenvectors just presented have a pathological behaviour at degenera- 
cies, since they become zero or linearly dependent and they do not form a basis. In the next we describe a procedure 
to renormalize them. 

3.3.1. Renormalized Alfven Eigenvectors 

We start by renormalizing the Alfven eigenvectors. From equation (|47p we see that the eigenvectors corresponding 
to the Alfven waves, ra,±, become zero for Type I and Type II degeneracies. To avoid this problem, let us first consider 
Type I degeneracy. In this case, Aa,± = Ae, leading to a = and B = 0, which is the reason for the Alfven eigenvectors 
to become zero. However, the quantity B/a is well defined for Aa,± ^ Ag, 



= tVs, (51) 

a,± 

and we can use this value to define the function B/a at Xa^± — Ag. Hence, we can obtain new Alfven eigenvectors just 
dividing the previous ones by a. These new eigenvectors are well defined even if we have a state in which the Type I 
degeneracy condition is fulfilled, and can be written as follows: 

ra,± = (e"0^^u'^(/.^6^ T^e'^p^.uf' <fh\ 0, 0)^. (52) 
Let us now consider degeneracy of Type II. The new Alfven eigenvectors just derived are proportional to 

= e'^p^.uf^'^h'. (53) 

This four- vector is orthogonal to u", 0" and 6". As discussed in the previous subsection, in the case of Type II 
degeneracy, can be written as 

= (^/.-f = (6,.f )^^^. (54) 
VCt -I- a"^ 

From this expression, it is apparent that 6" becomes a linear combination of and and, therefore, the four-vector 
e" is equal to zero for Type II degeneracy and so are the Alfven eigenvectors defined above. 

To avoid this problem, we proceed as follows. Using equation (I44p . we see that only the component of 6" tangential 
to the wave front, 6", has a nonzero contribution to the four-vector e", i.e., 

= e'^f.^.u^^-^hl (55) 

If we are considering an Alfven wave propagating in the x-direction, the vectors u", r" = (0,0,1,0) and r" — 
(0, 0, 0, 1) form a basis that we can use to decompose the vector 6": 

hi = gir; 53^" + 34^^"- (56) 

Substituting this decomposition of 6" in equation (|55p . we obtain 
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= gia'l + g2a!^, (57) 

where 



"i = 475"^ -^^^y = W{v^, Xv\0, 1 - Xv^, (58) 



«2 - e%5^^^^^! = -PF(t;^ At;^ 1 - A^;^0). (59) 
Coefficients 51 and 52 follow from ((57|) once is computed from its definition (equation [53|) . 

« = W + T^^') • (»' 

Vectors a^, 03 are clearly well defined through Type II degeneracy. On the other hand, since they form a basis of 
the orthogonal space to u" and z/" they can be used to express, as we have done with the vector e'', the tangential 
magnetic field bt as a linear combination of them, 

6^ = Ci^ + Caa^ (62) 
Coefficients Ci and C2 are related to gi and g2 as follows, 

Ci = ^i^i^±^Ty(l-Az;-), (63) 
aiia22 — 0^12 

C2 = -^i^il±^V^(l-Az;-), (64) 
aiia22 - 

where an = a'^ai^, ai2 = ol'^ol2^ and 022 = 0^2 0^2/^- From here it is easy to check that the necessary and sufficient 
condition for the tangential field component 5" to vanish is gi — g2 — 0. We use this fact to renormalize the Alfvcn 
right eigenvectors. If we divide by \/gf+g^, the right eigenvectors for Alfven waves read 

ra,± = (/i< + /2a^,TV£(/i< + /2a^),0,0)^, (65) 



where 



V 9i + 92 

and we take, at points where gi = g2 = 0, which corresponds to the Type II degeneracies, the following prescription 

These renormalized Alfven right eigenvectors are a linear combination of the ones proposed bv iKomissaro 3 (HMD 
for the Type II degeneracy case. However, contrary to the Komissarov's choice, our expressions are free of pathologies 
not only in the Type II degeneracy but also in the Type I degeneracy cas e. Note that the re normalization of the Alfven 
right eigenvectors just described is the relativistic generalization of the IBrio fc Wul (119881 ) method for classical MHD 
and the limiting valuse of /i^2 are chosen as in that work. 

3.3.2. Renormalized Magnetosonic Eigenvectors 

Let us consider the eigenvectors corresponding to magnetosonic waves. Again, if B^, = (Type I degeneracy), 
provided that a(A) and B{X) vanish for A = Ae, we have that As,± = Ae are solutions for the characteristic wavespeeds 
of slow magnetosonic waves. Since a = and B — 0, the corresponding right eigenvectors rs,± defined in (f48|) are 
zero. This pathological behaviour is removed if we divide the eigenvectors by and take the appropriate value of the 
quantity B/a when a — 0. A simple calculation establishes that for the magnetosonic eigenvalues Xm.± the following 
relation is fulfilled 




-M-^-Ag^ (68) 
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where m = s, /. 

After dividing by and taking into account equation (|68p . the new eigenvectors for the magnetosonic waves read 



C ' ph 



M 



where we have further divided by (c^ — 1)/Cs. 

These magnetosonic eigenvectors are weh behaved at Type I degeneracy as we can see by direct substitution of 
a = in the above equation, indicating that the pathological behaviour of the slow magnetosonic eigenvectors has 
been removed. Actually, for Type I degeneracy, the slow magnetosonic eigenvectors are given by 



which are a linear combination of the ones proposed bv IKomissaro ^ ([19991) . 

We now turn to Type II degeneracy. First, we note that the above magnetosonic eigenvectors, equation (|69p . have 
all their components proportional to \ht\ = y/bfhu' (see Appendix [B|) . Hence the components of the eigenvectors 
associated to characteristic waves for which — (0,0,0,0) (Type II degeneracy) will be equal to zero. To avoid this 
problem we can renormalize them by dividing by \bt\- If we denote by 

r™.± = (e^L^C,0)^ (71) 
the renormalized eigenvectors, from the expressions in Appendix [B|, we have 

e = u 2(n^ 2N +au)-[-] -7-77-1 , ("^^^ 



a2 - (G + a^)cl 



C = - .. ^r./AM . (74) 



Taking into account equation (p^ . the quantity b'^/\bt\ appearing in the definition of and L" can be written as 

bt _ (/iai2 + /2a22) ai - (/laii + /2ai2) 



[(aiia22 - q;^2)(/i"ii + 2/i/2Q!i2 + /|a22)] 



1/2- 



(75) 



where fi,2 are given by equation ((66|l with /i 2 = l/\/2 if gi = .92 = 0. 

The renormalization just described can be applied to the four magnetosonic eigenvectors: corresponding to two 
slow and two fast magnetosonic waves. However, in order to maintain the independence of the eigenvectors when 
we approach the degeneracy, we only apply the above renormalization to two magnetosonic eigenvectors, one of each 
class, those whose eigenvalues are closer to the Alfven eigenvalue of the same class. For them, we use expression (iTll) . 
For the particular subcase Type II', corresponding to Aa,± — ^s,± — ^f.±i ^a.+ 7^ -^a,-, for which the denominator of 
equation (j74p vanishes'^, we take C = and the problem is removed. 

For the remaining two magnetosonic eigenvectors we obtain renormalized expressions by dividing the eigenvectors 
by pha? /G — b . After some algebra, if we denote by Yrn.± = {e^ , L'^ ,C, 0)'^ the renormalized eigenvectors, we get: 

^ This result is easily derived from the following conditions: i) for Alfvc5n eigenvalues, a^/B^ = £; ii) for the degenerate eigenvalues in 
the Type II degeneracy, 3^/(0^ + G) = b^; iii) the definition of Type II' degeneracy, = b"^ /£. 
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ph 



Gb- 



G J pha^ - PG' 



(77) 



C = -1, 



(78) 



with the following prescription for the particular degeneracy subcase Type II' (for which a /G = b / ph) , 



0. 



pha? - b^G 

This procedure guaranties to have a complete set of right eigenvectors linearly independent for all possible states 

3.4. Renormalized left eigenvectors 



(79) 



We start with the vectors, 1,' 
Entropic wave: 



Alfven waves: 



le = (0„,0„,0,U° 



la,± 



b° ] ea,f3-,s^^un^ 



(«/a)a,± / 



\ 



where {B/a)a.± is given by equation ([5T|) . 
Magnetosonic waves: 



(80) 



(81) 



B 



b"\+bj- + 2a]{b^- 



m,± 



B 



B 



6" ~ -u'> ] b, 
a 



pha \ cj 



V 







a 



(82) 



where {B/a)m,± is given by equation ((68|l . 

Following the same procedure we have used to renormalize the right eigenvectors, we can obtain left eigenvectors 
well behaved for degenerate states. Alfven left eigenvectors take the form 

* It follows easily from i) — (G + ii^)c^ = 0, and ii) the definition of Type II' degeneragy, = b'^ /£. 

^ These vectors are related with those presented by Anile (1989) for the problem = 0, with (/)^ = {—A, 1, 0, 0), through 1 = 1.4*^. 
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\ 







(83) 



where fi,2 have been defined in equation (|66|) and we take again /i^2 = 1/^2 when 51 = (72 = 0. 

The magnetosonic left eigenvectors given by ([5^ are not weh defined for Type I degeneracy since a = B = 0. To get 
rid of this problem we multiply the eigenvectors by a so that 



A6^ 



m,± 



(/)^(-afe" + Su") + {acjP - Gu^)K 



(84) 



a0°-a2uO( - -1 ) + ^ 



5° 



V / 

4. A NOTE ON RELATIVISTIC MHD CONVEXITY 

In a convex system, all the characteristic fields are genuinely non-linear or linearly degenerate. Brio & Wu (1988) 
noted that the equations of classical MHD are non-convex since at degenerate states magnetosonic waves change from 
genuinely non-linear to linearly degenerate. We have checked the non-convex character of the relativistic MHD and 
found that the relativistic counterpart is also non-conv ex with the magnetosonic waves changing character also at 
degenerate states. We refer the reader to lAntonI (|2008[ ). where an analysis of the characteristic fields of RMHD in 
covariant variables is presented. 

The non-convex character of the classical (and relativistic) MHD equations is source of several pathological be- 
haviours, as the development of the so-called compound waves. The shock tube problem discussed in Sect. 18.1.21 
(proposed by Brio & Wu 1988 in classical MHD and adapted to relativistic MHD by van Putten 1993) displays the 
propagation of one of such compound waves. 

5. RIGHT EIGENVECTORS IN CONSERVED VARIABLES 
5.1. Transformation matrix 

The renormalized eigenvectors obtained in Sect. 3 are now transformed to conserved variables. To this aim, we must 
construct the transformation matrix between the set of covariant variables, U, and the set of the conserved ones, U, 
namely (i9U/i9U). At this point, let us recall that our aim is to build up a Ricmann solver based on the spectral 
decomposition of the flux vectors Jacobians of the system in conservation form. Since the Riemann solver is used to 
compute the numerical fluxes for the advance in time in a dimensional-splitting manner, only a onc-dimcnsional version 
of system ([T]) need to be considered. Consistently with the choice made in Sect. 3, we particularize our discussion to 
the x-direction. Along this direction, the evolution equation for is dB'^ /dt = 0, and it can be simply removed 
from the system. Hence, the desired spectral decomposition will be directly worked out for the shrunk, 7x7 Jacobian 
of the flux vector along the x-direction. Accordingly, in what follows the set of conserved variables will contain only 
seven variables and the aforementioned matrix will be of dimension 7 x 10. Its elements, corresponding to the partial 
derivatives of the conserved variables with respect to the covariant ones, are the following^ 

^ The transformation matrix is not unique since the constraints can be used in different ways to give different functional dependences. 
However, the resulting transformed eigenvectors are independent of the matrix used. 
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f p 

£uy £u° ~2b°u^uy-by 2b''u°uy 

£u'- £u° -26"m°m^ - 6^ 26^^°^^ 

2Eu^ -26"M"ii" - feO 26^(m°)2 - 6^ 

b?' -6° -wf 

V 6^ -&° -u^ 







26?^t 



26^mOu^ dp{ph)u^u^ ds{ph)u°u^ 



2byu°u^-b° 2b^u^u^ dp{ph)u"uy ds{ph)u°uy 





2byu°u'' 


26^uOm^ - 




2by{u^Y - 


by 2b'{u°f 




u° 












where 







,0\2 



,0\2 



ds 



5.2. Right eigenvectors in conserved variables 
The right eigenvectors in conserved variables, R, are computed as foUows 



R 



dV 



(85) 



(86) 



where r are the corresponding (renormahzcd) eigenvectors in covariant variables defined in equation (j46p -entropic 
vector-, equation (|65p -Alfven vectors- and equation (j7ip -magnetosonic vectors, with the corresponding definitions 
of e'^, L'^ and C, for each pair of vectors-. 



5.2.1. Entropic eigenvector 
The entropic eigenvector in conserved variables is 

Re = u°{d,p, w'dsiph), uyOsiph), u^a,(p/i), u^d^iph), 0, 

5.2.2. Alfven etgenvectors 



(87) 



Alfven eigenvectors are 



where 



R 



a,± 



/iVa, 



1,± 



2,±, 
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Va,l,± = 



Va,2.±=- fw^'u^i V^^^^U^ 

V ± VSuyu'^ 



(90) 



Vectors Va^i^± and Va,2.± are well defined in Type I degeneracy. In tlie case of tlie Type II degeneracy, coefficients 
/i and /2 in ([55]) must be taken equal to 1/ \/2. 

5.2.3. Magnetosonic eigenvectors 

Finally, in terms of e", {v — 0,x,y,z) and C, the components of the magnetosonic eigenvectors in covariant 
variables defined in equations (|72p - (j74p respectively, the magnetosonic eigenvectors in conserved variables are given by 

/ pe° + Cu^dpp \ 

e(u^e^ + u°e^) + 2u^u'=baL°' - b^'L" - b^L"" + u°u''Cdp{ph) 
8{uye° + u°ey) + 2u°uybcL°' ~ byL° - b°Ly + u^uyCdpiph) 
£:(u^e° + + 2u°u'-baL°' - V^L^ - b°L^ + u°u^Cdp{ph) 

2£u^e^ + {2{u^Y - l)6„L" - 26"^° + {dp{ph){u^ f - l)c 

yye^ - b°ey - uVL" + u°Ly 



R 



(91) 



According to the renormalization recipes given in Sect. 3, the components of the magnetosonic eigenvectors corre- 
sponding to the eigenvalue closer to the Alfven one within each class, are 



R 



albt 



m,±, D 



/i(a2 -c2(G + a2))' 



^0 + „,0)_M^^ 



1^*1 ^.^t^'^l^y dpp^ 



a2 - (G + a2)c2 



(92) 
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R 



G J ia2„c2(G + a2) 
a 



-uVC9p(p/i), (93) 



Rm.,±.,Si — I 1 



G J ia2„c2(G + a2) 



[a^'il-cl) + 2u'clG) 



+ nrr + irr + irr r + u^u^cdpiph), (94) 



a2 



(95) 



G |6t| \bt\ 

where in the preceding expressions, j stands for y and z, and a, G and (/)^ are computed for the corresponding 
magnetosonic eigenvalue, \m,±- 

These eigenvectors are aheady well defined in Type I degeneracy. Quantity b^ /\bt\, depending on /i^2, was defined 
in equation ([75]) . In the case of Type II degeneracy, coefficients /i and /2 must be taken equal to 1/ v^. 

For the components of the two remaining magnetosonic eigenvectors we have 

7? ° , Q\ (g/a)m,± b^tG Oq ,c,n\ 

^-■^■ ^ = hiG + a^)c^^ '^ h pha^ -h-G -"" ^^^^ 

Rm,±,s^ = ^(a(l-c2)(u-0O + uV^) + 2uVc2G) 



p/ia2 - 62G 



6^6° + 6^5^ - (6^u" + 6°u^) ( - j ) - u%^c)p(p/i), (98) 



'U? f 

^m,±,5. =^(a'^"(l-c?) + 2u"c2G 



Gc2 



fe^foO + foOfoJ' - (b^u" + 6?M^) j - u''u^9p(p/i), (99) 



p/ia2 - 62(7 

. = ^ ( 27."(a0°(l - cl) + u°c2g) + - c2(G + a") 



R 



'771, lb, r — 



26? (G 



p/ia2 - 62G \ Vo / m,± 



(100) 



1-^ ((u-A™,± - ^.°)6^" + Gu=bi) , (101) 

where, again, j = y,z, and a, G and </)^ are computed for the corresponding magnetosonic eigenvalue, A„i,±. In the 
special subcase of Type II' degeneracy, 6J' l{pha^ — b'^G) must be taken equal to 0. 
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6. LEFT EIGENVECTORS IN CONSERVED VARIABLES 

6.1. Transformation matrices 

The calculation of the left eigenvectors in conserved variables is far more involved than that for the right eigenvectors. 
One possibility would be to obtain them by direct inversion of the matrix of right eigenvectors. However, in order to 
obtain tractable expressions, we have proceeded in a two-step process starting from the left eigenvectors in covariant 
variables presented in Section 3. First we get the left eigenvectors, 1, in the so-called reduced system of covariant 
variables, V = {u^,uy,u^,by,b^,p,p), through the transformation 



1 = 1 




(102) 



where (9U/9V) is the 10 x 7 matrix built with the partial derivatives of the covariant variables as functions of the 
covariant variables in the reduced system. 

In a second step, the left eigenvectors in the conserved variables, L, are obtained from those in the reduced system 
of covariant variables through 



(103) 




L = 1 

We give now the expressions of the two matrices. For the first one we have 

/ du° du° du° 




du^ duy du^ 







1 























1 























1 














db° 




db° db° db^ 










duy du^ 


dby db'^ 


d¥ 


db"" 




db"" 












duy du^- 


dby db^ 











1 























1 























1 







V 














(1). 


(ds 



(104) 



In the previous expression. 



= -fT i = x,y,z, 



db° 
W 



(u0)2 - (u-)2 



db° 
du^ 



(105) 
(106) 
(107) 



db^ _ 1 



y,z, 



(108) 
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The second matrix, {dV/dU), involves the derivatives of covariant variables with respect to the conserved ones. 
In order to obtain general expressions, we use the variable Z = ph{vP)'^, introduced in Sect. 2. Once the function 
h = h{p,p) is provided through the equation of state, the full system is closed and the partial derivatives of Z can be 
calculated. 

As a first step, we write vP, 6° and as functions of the conserved variables and Z: 

= -J., (112) 

{Z + B2)2 - S2 - ^ '^J {2Z + B2) 

b' = {S-B)^, (113) 

'=^+ Z + B2 (114) 

From these expressions, the corresponding partial derivatives with respect to the conserved variables can be written 
in terms of the derivatives of Z. Then, the covariant variables in V are written in terms of the conserved variables, 
W and 6^, and Z 

u = z + B"^ {i = x,y,z), (115) 
V = ^-^^ {i = y,z), (116) 

, = ^. (118) 

The elements of the transformation matrix are obtained by derivation of the previous expressions and substitution, 
leading to: 

du' -vP ( n , b%^\ dZ 



u 



« / ^ b%'\ dZ f i B^b \. , , 

5;._UV + — — + , (120) 



du' -u" / , b%'\ dZ 

« " + ^ Ur. (121) 



dr Z + B2 V Z J dr 
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dD Z + B2 ' 



(123) 



dS3 Z + B'^y\ V ' I Z ) dSi 

[b, (i + ^^) +50^^.^ _BS + ^<5i|, (124) 



+ 



^] + ( 1 + ^ ) + ^ (50.' - B^) - 2i., (."5'^ - ^) I , (126) 

z + az 



dD Z + 9Z) 

9p _ Z + 1 

z + 62 



(127) 

B^uj - ) , (128) 



dr Z + B2 ar 



- 1, (129) 



dp 1 j , 2n 5^ ^"-S-j 52^^. ;l 



" ^-(.y-mi, (131) 



az? 1*0 mO(Z + B2) V Z J dD 



aS-J u"{Z + B^)\\ ' ' Z J dSj V ^ 



5t uO(Z + B2) V ' ' Z ; c)r ' 



l-{uy-^ i^ + 2BAl-iur +b^>r^+uA . (134) 



mO(Z + B2)1 V ^ ' Z J dB3 '\ ^ ' J \ Z 

6.2. Left eigenvectors in the reduced system of covariant variables 
6.2.1. Entropic eigenvector 

Multiplication of the entropic eigenvector in covariant variables and the transformation matrix ((9U/9V), defined 
in (|104p . leads to the entropic vector in the reduced system of covariant variables 

le- fo,0,0,0,0,«°('^) ,u'(?^) V (135) 



dpjp \dpyp. 
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6.2.2. Alfven eigenvectors 

Performing the same operation with the renormalized Alfven eigenvectors in covariant variables, we get 



where 



and 



Wa,l,± = 



la,± = /lWaa,± + /2Wa,2,±, 

£u"uy ± V£{hyu° - 2b°uy)juyu° 

u''{±y^uy + by) 
tV^(i + K) 

V J 



( 2{-£a + b°B/u")uy 

£:(i + (u^)2) ±V£(^b°(u" ~2{uy)^/u°^ +byuy -b^'u'' 

- Eu^u'^ ± V£{b'u" - 2b^u^)^uy/vP 
-¥-u'- =F fl + {u'^f 



V 



) 



(136) 



(137) 



(138) 



Vectors Wq_i.± and Wa^2.± are well defined in Type I degeneracy. In the case of the Type II degeneracy, coefficients /i 
and /2 in (|136p must be taken equal to l/\/2- 

6.2.3. Magnetosonic eigenvectors 

Finally, multiplication of the left eigenvectors in covariant variables and the transformation matrix (9U/9V), defined 
in (|104p , and further algebraic manipulation lead to the following expressions for the components of the magnetosonic 
eigenvectors in terms of A, gi, (?2 and bf (see Anton 2008), 



L,±,„^ = ^ {{uy - inn') + 2BX{G + a^)) , 



(139) 



u 



/2S^A,„,±6°(G + a2) A 



g,UG + 2a')\b^>- 



B 



u 

a 



i^x 



6-A„,,± + 6° L 



(140) 



m.zb 
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- _ f 2B^Xrn,±b'}{G + a^) A , 

g2UG + 2a')(b°^(-) uA - Xm,± + bA , (141) 



a 



-5i("°-A™,±u"), (142) 
lm,±M = -52(u° - A™,±m"^), (143) 

L,±.p - ^ v(^"-±"" - - (-) ^t, (144) 

L,±,p = 0. (145) 

The inspection of the previous components allows us to conclude that all of them are proportional either to A, gi , 
32 or which are all zero in the Type II degeneracy. Then the magnetosonic left eigenvectors just obtained are not 
well suited for Type II degeneracy and must be renormalized. 

In a first step, we eliminate A taking into acount that, for the magnetosonic eigenvectors^ 



a\G + a^){l~cl) 
a2-(G + a2)c2 



^ = -^^7T^^7IT:#^^t • (146) 



Similarly to the case of the right magnetosonic eigenvectors, those eigenvectors, one of each class, whose eigenvalues 
are closer of the Alfven eigenvalues, are then divided by \bt\. The new eigenvectors have components: 

- ^ (^^^i^(<"")^ - <"-)0 - . <-) 

- uy{G + a^) ( {l~cl)\b,\au- 6° 

- u^{u^~\„,^±u^) \a^~WTa^}^ -^"'^ -u)-2B A™,± — 

+ |L / (G + 2a2) {b^-(-\ uA - 6-A„.,± + 6°) , (148) 



6-A™,±+5n, (149) 



lm,±,b« = -^{U° - \m,±un, (150) 

r„,±,b. = -T7^(w° - A™,±u"), (151) 

''' The way to obtain this expression is a bit tedious: substitute in equation IIB3t . valid for the magnetosonic eigenvectors, by 
= (G + a?)({ph + b^)a? — A)/G. Previously, obtain this expression starting from the definition of , equation II40I I. then substitute 
in £ in terms of its normal and tangential parts according to equation I IB2I I. and work out the value of B^ from the resulting expression. 
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lm,±,p = 0. (153) 

In the previous expressions, follows from equation (|75p . and 

,1/2 

(154) 



/l.2 ((^i" " Xm,±un' - (1 - A,2„,±)((U^)2 + (U^f)' 



1^*1 (/|«22 + 2/1/2^12 + /2an)'^' 

In the Type II degeneracy, coefhcients /i and /2, in ([75)1 and p54p . must be taken equal to l/\/2- 

For the remaining two eigenvectors, those of each class whose eigenvalues are farther from the corresponding Alfven 
eigenvalue, we propose to divide the original components by the expression 6(/(a^ — (G + 0^)0^), and make the 
substitution 



a2-(G + a2)c2 G(G + a2)c2' 
valid for the magnetosonic eigenvectors. The rcnormalized eigenvectors are 



(155) 



W^l(<i^((.T-K.').^), 

uy ( {l-cl)au% 2S-A„,±60G^ 

I 2 ^'^ ^m,± — U j 



pha^ - b^G 



(G + 2a2) U"- (^^^ uA-b^X^,± + bA, (157) 



- ^ (l~c2)a^,^ ^ 2i3-A™,±60G 

'in,±,ti» — r, / n ^ TT 5 I," '^m,± — " J 



uO(uO- A™,±M-) V c2 ' ' pha^-b^G 

+ ;i™{<= + ^°'> ('"-(!).„,, "°) + (158) 

L,±,b« = ■ 2^ 1.9^ - Am,±M''), (159) 

p/ia^ — o^G 

p/ia^ — o^G 

- _ (A^,±7.- - ^°)G G2 6° 

pMG + a2)c2 p/ip/ia2-62G ^^^^^ 

L,±,p = 0. (162) 

In the case of Type II' degeneracy, quantity b^ /{pha^ — b^G) must be taken equal to 0. 

We note that in all the expressions of this section, quantities depending on the eigenvalues, have been computed in 
terms of the corresponding Xm,±- 

6.3. Left eigenvectors in conserved variables 
6.3.1. Entropic eigenvector 

The last step involves the multiplication of the entropic eigenvector in the reduced system of covariant variables by 
the matrix (9V/9U), whose elements have been defined through the expressions (|119p - (|134l) . The resulting eigenvector 
in conserved variables is 
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^ 1 9s 1 



Wdp Z + B^ 



with i = x,y, z, 

- " - V (2?: . .)) I - - 1^) I)' 

with i = y,z. 

6.3.2. Alfven eigenvectors 

The same operation performed now with the Alfven eigenvectors leads to the corresponding eigenvectors in conserved 
variables 

L„,± = /iW„,i,± + /2W„,2,±, (167) 
where the components of vectors Wa,i,± and Wo,2,± are 

dZ 

Wa,l,±,D = Ca,l,±Q^, (168) 

BZ 

W^a,i,±,s« = Cax±^ + (B^u^ - 6°^-)^ - u-hy{Byu' - BV) 
+ {£vP ± \/f 6°)u^(m^ - 2a) ± ^{u'^B'' - u^B'' - uyu'^iu^By - uW)^ , (169) 

T^a,i,±,s» = Ca,i ±^ + {B\y - iPBy)— + h%yu' - uybyiByu' - b^u^) 
Bt>y u" 

± VS (^b^u^u^ + (u^)2) {uyB' - u'By)^ , (170) 
BZ 7/^ 

Wa,i,±,s' = Ca,i,±^ + {B^u^ - b°B')-^ - b%yuy - u%y{Byu' - B^uy) 

+ £{uy{u° - Aa±u^) ± V£(b°u''a + b'u'u° - u'uy{u'By - uyB')), (171) 

BZ 

Wa,l,±,r = Ca,l,±^-u\Z + B^), (172) 

_ BZ \ B^hy-b^sy / 1 \ 



+6^(B^.^ - Byun + 25^^ ^) 



+ ± b°V£) + (1 - au'=){h%y - 2{u^fBy) J - ^6^w-a) 

+by{B'uy - Byu') f ^ + 2^^ ^'^^' ^ 

+(6:^0 ± ) + (1 - au')(fu' - 2{u°fB') J - ^h^u'c^ 

- 2iu^fB^) ^ - (i.^ + ^B^uy - Byu^uy) + 2^^!^) } . 

Wa,2,±,D = Ca,2,±-Q^, 

8Z yy 
+ {£u° ± V^6°)w^(w^ - 2a) ± - u^B^ - - u^B^)) , 
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+{b\y-2iurBy)^^^ 



(173) 



(174) 
(175) 

(176) 

(177) 



'a,2,±,Sv = Ca,2,±^ + (B^U^ - 6°^^)^ - b°b'u' - uyb^B^uy - By u') 

+ £{u°f{u° - Aa±u^) ± \/f (6°u^a + fo^u^u^ - u''uy{uyB'' - u'^By)^ , 

Wa,2,±,s^ = a,2,±|^ + {B\' - b''B')% + b%'uy - u'b'{B'uy - Byu') 

± V£{b'uyv° + + ^'^'^^) ("^^^ " "^^^)) ' (1 ''8) 

dZ 

Wa,2,±,r = Ca,2,±^-uy{Z + B^), (179) 

w.,2,.,s. = a,,,— + + (2 - p^j zBy] uy 

+b'{Byu' - B'uy) + 2B^ -^~y ^ 

+ (f^o ± 6°V5) + (1 - au-) (b\y - 2{uyBy) ^ - 

T {^^^^ (l + + - ^) {{ur - K^u'u^ - 



(180) 
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+ (£7/° ± \P^f£) ( + (1 - au-) ( b^/' - 2{u'> fB') ^ - ^b^'uya] 
\ u" \ / m" Z J 



+ (5^1/ - — - + 2B^ r • (181) 



In the previous expressions, 



(182) 



Ca.i,± = -Su°u'{l - aw^) - {8u° ± V£b")—{b''{l + aw^) + u^ib^uy - byu') ~ 2b'=au'j 

+{z + b^)u'' - by{v-uy - byu") [{uy - 1) 

Ca,2,± = ~£u°uy(\ - au"") - {£u° ± VSb")—{by{l + au"") + u'{byu' - ¥uy) - 2b''auyj 

+ {Z + b^)uy - b'ibyu' - 5"u^)((u°)2 - 1) 
±V£ ju^' (^{u'Yu^B + ^"^^^" ^ - ((w°)2 - 1) (b?' + u^byu' - b'uy)^ I . (183) 

Vectors Wa.i,± and Wa_2.± are well defined in Type I degeneracy. In the case ol the Type 11 degeneracy, coefiicients 
/i and /2 in p67p must be taken equal to 1/ \/2. 

6.3.3. Magnetosomc etgenvectors 

We present now the magnetosonic eigenvectors in conserved variables from the transformation of the corresponding 
eigenvectors in the reduced system of covariant variables before the renormalization process for Type 11 degeneracy, 
those defined through expressions (|139p - p45p . We start by defining the quantities 

H„,,±^{G + 2a')lb'>-(-] u° ) -(5^A^,±-0, (184) 



a2 - (G + a2)c2 

„2 i\i-r< I „2\„„,0 „.x , " ^ 111 f„.0\2 I \ /\ „.x „,0\ 



X |(cj-l)(G + a^)au" i^u^ + ^ j + |^(u")^ + — j (A„,±w- - u")G 
,0 [2(G + a2)B-uO / ^ogx h''\^fB 



52 i w 



Now, the components of the magnetosonic eigenvectors read, 



(185) 



Lm,±,D = Cm,±-gJ^, (186) 
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+ 52|«°//m,± (u^u' + ^) - - \m,±u' ) [ii' [h^'vP + B-^ - B^u^y^, (187) 



_ 9Z 

+ 5i|u°F„,± (^1 + {uyf + - {u° - Xm,±u^)(uyBy^ + 6o(i + K)^)) | 

f/1 2N 2n5"-B^ B^u'-b°B\^ ^ p,^] 

+ 5i|w°ifm,± (^u'uy + - {u° - x^,±u^) (^uy (^b'u° + B'^ - Byu')^ I 

+ 52|n0if„,± (1 + {u^f + ^) - («o - Xm,±u^) (n'B^^ + fe" (l + {u^f) ) | , 



+ 

+ 6' 



X I (1 - cl)a{G + a2) (^^^ - 2B^w-wO^ 



(188) 



(189) 
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+syuy + + ^{b"uy - By)] - ( u\yb" + 

\ ph J J \ 



,0^2 

-By 

g2{H^^±{{b%y-2{uyByy + ^^-{u°-Xm^±un(^syu^ (^1 + ^ 

^(6V - B')^ - 2By (^uV6" + 1— Ij^i?-^ |, (191) 

r -r ' 



x{{l- c2)a(G + a") - 2B'u^u" 

I 2g-(G + a^) A„^^, ^ 5^ _ 2."i?^(A„,± + u%) 
^605^-526- ^^^^ ^^^,\fB 



ph u'^ 

-.2 



2(u") ) ( - ) G 



+ 91 i^H„,,±(^{b%^ - 2iuyB^) uy + - (^.0 - A™,±zi^) (^5^^.'^ (^1 + ^ 

+ ^(6°u^ - S^)^ - 2B^ + i^ij^B^^ I 



92 (Hr 



Our proposal of renormalization for Type II degeneracy is the following. For those magnetosonic vectors, one of 
each class, whose eigenvalue is closer to the corresponding Alfvcn eigenvalue, we divide by \bt\. Hence the following 
replacements in the components of the eigenvector (defined previously) must be performed 

^ (193) 



a2-(G + a2)c2 a2-(G + a2)c2' 

where the expressions that must be used in the case of Type II degeneracy are defined in ([75|) and ()154p . 

For the remaining two eigenvectors, those of each class whose eigenvalues are farther from the corresponding Alfven 
eigenvalue, we propose to divide the original components by the expression b1/{a? — (G + a?')c^g), and make the 
substitution 



62 _ pha'^-b'^G 

a2-(G + a2)c2 " G(G + a2)c2- 

Hence the following replacements in the vector components shall be performed 



(195) 



a2-(G + a2)c2 ' * pha^-b^G ' 
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^'"^ pha^-b^G ' pha^-b-^G ' ^ ^ 

where the expressions that must be used in the case of Type II degeneracy are defined in ([75)1 and (|154p . 

6.3.4. Normalization factors 

Finally, for the implementation in the numerical code the left eigenvectors just defined must be multiplied by 
normalization factors in order to fulfill the condition Rj — Sij . 

7. HRSC SCHEME FOR RELATIVISTIC MHD 

As established in the Introduction, HRSC schemes have become a standard tool in numerical relativistic hydrody- 
namics and magnetohydrodynamics. In this Section we describe the ingredients of our numerical algorithm, the most 
important being the linearized Riemann solver based on the renormalized spectral decomposition of the RMHD Jaco- 
bians described in the preceding sections. Although the specific expressions for the spectral decomposition have been 
worked out explicitly for the x-direction, symmetry arguments allows one to obtain the equivalent expressions for the 
remaining directions easily. The starting point are the equations of RMHD written as a first-order, fiux-conservative, 
hyperbolic system of partial differential equations ([7|). All the ingredients described in this section have been imple- 
mented in the code RGENESIS (Leismann et al. 2005), that we have used to perform all the numerical tests and 
applications presented in this paper. 

7.1. Integral form of the RMHD equations 

To apply this technique to the system ([7]) we need first to obtain an integral form of the equations. Integrating the 
system on the three-dimensional volume bounded by the surfaces S^,, 'Sx+Axj Sj^, Sj^_|_Ay, S^, S^+azi that connect two 
temporal slices, we obtain 
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I Y^dydz- I Y^dydz \ - I j F^dxdz- f F^dxdz) 

F^dxdy- j F^dxdy] , (198) 



where 



and 



1 



x+Ax ny-\-Ay pz-\-Az 



F° = — / / / F" dxdydz (199) 



AV = AxAyAz. (200) 

The carets appearing on the fluxes denote that these fluxes, which are calculated at cell interfaces where the flow 
conditions can be discontinuous, are obtained by solving Riemann problems between the corresponding numerical 
cells. These numerical fluxes are discussed below. 



7.2. Induction equation and divergence-free condition 

Independently of the particular Riemann solver we use for calculating the numerical fluxes at cell interfaces, the main 
advantage of the numerical procedure just described to advance in time the system of equations, is that those variables 
which obey a conservation law are, by construction, conserved during the evolution as long as the balance between the 
fluxes at the boundaries of the computational domain is zero. Although this is one of the most important properties a 
hydrodynamical code should fulfill, there is an additional condition we should demand to any magnetohydrodynamical 
code: at any time in the numerical evolution, the divergence of the magnetic field must be zero. This condition is 
not satisfied by construction if we use equation ()198|) to evolve the components of the magneti c field. In o rder to 
preserve the divergence free conditio n, we use the con straint transport method designed by Evan s~fc HawlevI (p988), 
and extended to HRSC schemes by iRvu et aTl (|1995D . The essential of this method is to use Stokes theorem after 
integration of the induction equation on the surfaces that bound our numerical cells. This allows one to obtain 
evolution equations for the magnetic flux at each cell interface in terms of the electromotive force around the contour 
defined by the boundary of the cell. Since the sign of the electromotive force changes if the direction of the contour is 
changed, the magnetic flux through the faces of a numerical cell will be constant in time. If initially this flux is zero, 
this condition will be fulfilled during the evolution and the divergence free condition numerically satisfied. 

Taking, for example, the magnetic flux ^s^, through the surface E2, defined by z = const., and the remaining two 
coordinates spanning the intervals from a; to a; + Ax, and from y to y + Ay, we have 
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= / B ■ dE. (201) 



The discretizcd form of the induction equation based on the method just described is the foUowing 

nidx\ (202) 



dt 



rii = EijkV^ B^ is the electric field with opposite sign, and the caret denotes again that quantities Vli are calculated at the 
edges of the numerical cells, where they can be discontinuous. At each edge, as we will describe below, these quantities 
are calculated using the solution of four Riemann problems between the corresponding faces whose intersection defines 
the edge. However, irrespective of the expression we use for calculating 17^, the method to advance the magnetic fluxes 
at the faces of the numerical cells satisfies, by construction, the divergence constraint. 

7.3. Spatial and temporal order of accuracy 

In order to increase the spatial accuracy of the numerical solution, the primitive variables defined in equation p7p 
can be reconstructed at the cell interfaces before the actual computation of the numerical fluxes. We use standard 
second order (MINMOD) and third order (PPM) reconstruction procedures to compute the values of p, p, Vi and B^ 
(i = 1,2,3) at both sides of each numerical interface. However, when computing the numerical fluxes along a certain 
direction, we do not allow for discontinuities in the magnetic field component along that direction. Furthermore, the 
equations in integral form a re advanced in time using the method of lines in conjunction with a second/third order, 
TVD Rungc-Kutta method (jShu fc OsheilfTOSl . 

7.4. Numerical fluxes 

As discussed in the Introduction, there are several strategies to calculate the numerical fluxes in HRSC schemes. 
Our procedure is based on the construction of a linearized Riemann solver after the original work of Roe (1981). The 
main idea of Roe's approach is to linearize the hyperbolic system of equations and to use the analytical solution of 
this linearized system at each interface to obtain the fluxes across the interfaces. These fluxes are then employed 
in the discretized equations to advance in time the variables. Linearizing the system amounts to choose a constant 
Jacobian matrix at each interface. In the original Roe's solver (Roc 1981), this constant Jacobian matrix is built to 
give the exact solution of the nonlinear Riemann problem if a discontinuous wave is located at the interface. However, 
it is possible to relax the conditions imposed to the Roe's matrix and to build simply the Jacobian matrix for some 
intermediate state between the two states separated by the interface, and still producing accurate solutions. Riemann 
solvers following this approach are commonly known as Roe-type Riemann solvers. In our numerical implementation, 
we use the arithmetic mean between the variables at each side of the interface to obtain the intermediate state. The 
variables we use to obtain the intermediate state are p, p, ^v^ , and B^ . Knowing the intermediate state we 
obtain the Jacobian matrix, and the eigenvalues and right and left eigenvectors in conserved variables as described in 
Sects. [3l [5] and [6] Contrary to the HLL strategy and its sequels. Roe-type Riemann solvers calculate the numerical 
fluxes consistently with the breakout of the original discontinuity in the full set of characteristic waves (full wave 
decomposition Riemann solver; FWD, in the next). Explicitly, numerical fluxes are computed according to 



F 



FWD 



1 



2 

where 



r(Ui) + F(Ufl) - I A^P) |5(f)R(p) 



(203) 



a 



ip) 



L(f)(UH-Ui), (204) 

\^P\ R(p) and L*^^^ being respectively the eigenvalues and right and left eigenvectors of the Jacobian matrix for the 
intermediate state. 

One of the problems of Roe-type Riemann solvers is that the entropy condition is not satisfied through transonic 
rarefactions. To solve this problem, an extra viscosity term is required at these points (Harten & H yman 1983; Ygi 

I1987D and this can be implemented by substituting the eigenvalues A'^''^ in equation (|203p by max (0,|A(p)|,(A(f) - 
Al^): (A^^) -A(f))y where A^^^ and A^^^ are the eigenvalues associated to the left and right states, respectively. 



In Sect. [3] we obtained the characteristic equation whose zeroes are the eigenvalues. In the case of the entropic and 
Alfven eigenvalues, the solution of the characteristic equation leads to analytical expressions which can be used directly 
to compute them. Magnetosonic eigenvalues, however appear as roots of a quartic equation which is solved numerically. 
To solve the quartic equation, we use a Newton-Raphson algorithm to obtain the eigenvalues corresponding to the fast 
magnetosonic waves using the light speed (-f 1 and —1, in the units of the numerical code) as seeds. Once obtained these 
two eigenvalues, the quartic equation is reduced to a quadratic one which is solved analytically. This procedure is very 
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accurate in the determination of the eigenvalues. Only when either the kinetic or the magnetic energy of the physical 
state is much larger than thermal energy, we find eigenvalues with large uncertainties (0.1 % or larger) satisfying the 
quartic within the machine roundoff error (10~^^, in double precision). The analytical solution of the quartic does not 
represent an improvement, since the large number of operations involved tends to increase the numerical error. 

The flux formula (|203p is used to advance the hydrodynamic variables according to equation (|198p and to calculate 
the quantities fi^ needed to advance in time the magnetic fluxes following equation (|202p . At each edge of the numerical 
cell, rii is written as an average of the numerical fluxes calculated at the interfaces between the faces whose intersection 
deflne the edge. Let us consider, for illustrative purposes, ilx- If the indices (j, k, I) denote the center of a numerical 
cell, an a;— edge is defined by the indices (j, k + 1/2, 1 + 1/2). By definition, H.^ — (w^i?^ — w^B^). Since 

F^iB^) ^vyB" -v^Sy (205) 

and 

F'{By) = v'By -v^B', (206) 
we can express in terms of these fluxes as follows 

^xj,k+l/2,l+l/2 = ■^i^lk+l/2d + ^j.k+l/2d + l ~ ^lkd + l/2 - Pj.k+l.l+l/2]> (207) 

where Fy{F^) refers to the numerical flux in the y {z) direction corresponding to the equation for B^ {By). 

In our numerical scheme, we need also to know the value of the magnetic fleld at the center of the cells in order 
to obtain the primitive variables after each time step (see below) and to compute again the numerical fluxes of the 
other conserved variables for the next time step. If -BJ±i/2 k i x-component of the magnetic field at the interface 

(j ± 1/2, k, I), then the x-component of the magnetic field at the center of the (j, k, I) cell, BJj, is obtained by taking 
the arithmetic average of the corresponding fluxes, i.e. 

Bj.k,l = 2(-^j-l/2,fc,i^'S'j-l/2,fc,i + Bj + i/2,k,l^^j+l/2,k,l)/^^3.k,h (208) 

where AS'Jj_j^^2 k i area of the interface surface between two adjacent cells, located at Xj±i/2 and bounded 

between [yfc-1/2, yfc+1/2] and [zi_i/2, zi+i/2]- Analogous expressions for ^yj+i/2,k,i+if2 and ^zj+i/2,k+i/2,i, and B^jj^ ^ 
and Bj f. ; can be easily derived. 

7.5. Recovery of the primitive variables 

After every integration time step, primitive variables are recovered following the procedure described in Sect. 2, as 
in Leismann et al. (2005). 

8. NUMERICAL TESTS 

In this Section, we evaluate the performance of the Riemann solver presented in Sect. 17.41 bv means of a number 
of one dimensional shock tube problems, as well as, some two dimensional standard tests previously computed in the 
literature. 

8.1. One dimensional tests 

For easier comparison with previous works in the field, we perform the same one dimensional tests as proposed by 
MUB09, with initial states to the left and to the right of the discontinuity given in Table 1. The initial jump is set 
at a; = 0.5 in a domain spanning the range [0, 1], which is covered by uniform numerical zones (see Table 1). A 
constant 7-law with 7 = 5/3 is assumed for the equation of state. We perform all the tests with a CFL number equal 
to 0.8. By default, RMHD equations are integrated using a second order Runge-Kutta integration with no spatial 
reconstruction in order to assess the performance of our FWD Riemann solver. We will refer to this basic setting as 
our first-order scheme (in analogy with MUB09). 

To show the num erical accuracy of our Roc- type Riemann solver, we also compute the numerical solution with both 
HLL ( Leismann et al . 2005) and HLLC (Mign one fc B odo solvers. Tests with the newer HLLD solver are not 

performed but our results are directly compared with those presented in MUB09 for the same numerical tests. The 
discrete errors in L-1 norm are computed as 

2—1 

where Vi is the first-order numerical solution (density or magnetic field) of our scheme and Vf''^ is the analytic solution 
at Xi, computed with the exact Riemann solver of .Giacomazzo Sz Rezzolla (2006.) . 
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8.1.1. Contact and Alfven discontinuities 

In an isolated contact wave, only the rest-mass density is discontinuous. In Fig. [1] (left), we show the solution of 
the first Riemann problem of Table 1 at i = 1 as computed with the HLL, HLLC and FWD solvers. Both HLLC 
and FWD solvers show no smearing of the initial profile, while HLL displays the largest numerical difussion. This 
behavior is expected for the HLLC solver, which is specifically built to include a contact wave in the numerical solution 
(jMignone fc Bodoll2006[ ). Also in the FWD solver the use of the full spectral information of the RMHD system yields 
perfectly sharp contact waves. 

We set up an isolated rotational wave as in MUBG9, and show the variation of By across it. The right panel of Fig.[l] 
shows that both HLL and HLLC solvers are unable to capture the initially specified jump (see Table 1), while the 
FWD solver captures it with very little diffusion. In this regard the Roe-type solver compares reasonably well with 
the HLLD solver of MUB09, which exhibits no diffusion at all at this kind of discontinuities. 

We shall note that, in spite of the fact that HLL smears both contact and rotational discontinuities, in practical 
applications, the numerical diffusion is not as bad as one may guess by extrapolating the results of Fig. [T] If we use 
a third order accurate scheme both in space and time (using PPM for the spatial interpolation, and a third order 
Runge-Kutta time integrator ^, the numerical diffusion is drastically reduced (Fig. [21). Even HLL is able to resolve the 
contact discontinuity in a couple of points (Fig. [2] left), and both HLL and HLLC need only 4 to 5 numerical cells to 
resolve the rotational wave. 

8.1.2. Shock Tube 1: The Brio & Wu test 



iBalsaral ()2001[ ) proposed a relativistic extension of the classical magnetized shock tube of iBrio fc Wul pSl) 



which has also been considered in a numb er of recent papers (e.g., iDel Zanna et al.l 120031 : iLeismann et al 



I2OO5I : 



iMignone fc Bodol 120061 : iMignone et al.l [2009? ) . In this case, the adiabatic index is 7 = 2, and the initial piecewise 
discontinuous state breaks into a left-going fast rarefaction, a left-going compound wave (see Sect. |4|), a contact 
discontinuity, a right-going slow shock and a right-going fast rarefaction wave. 

The final state of this test using the first order scheme is shown in Figs. [3] and 21 Our FWD solver performs clearly 
better than HLL or HLLC, and the quality of the results compares fairly well with the HLLD solver of MUB09. The 
superior quality of the results using the Roe- type solver is more obvious in the central part of the Riemann fan (Fig.lU, 
which is much better resolved than with the HLL or HLLC solvers. A more precise quantification of how better is 
the overall solution can be seen from Fig. [5l (upper left), where the Ll-norm errors associated to the FWD solver are, 
systematically ~ 54% and ~ 40% smaller than those associated to the HLL and HLLC solvers, respectively. 

The computational cost of the FWD solver is higher than that of the other two solvers, being the ratio of CPU times 
for this test and the first order scheme ^hll : ^hllc • ^fwd = 1 : 1-02 : 2.92. Using a the third-order scheme, the CPU 
times scale as ^hll : ^hllc '■ ^fwd = 1 : 1-05 : 2.48. 

8.1.3. Shock Tube 2: Non-planar Riemann problem 

I n this test, proposed bylBalsaral (|200lD and also performed by, e.g.. ILeismann et all (|2005[ ): [Mignone fc Bodol (|2006f) . 
and IMignone et all (2009*), the transverse magnetic field rotates across the discontinuity by ~ 0.557r, and the Riemann 
fan is composed of three left-going waves (fast shock, Alfven wave, and slow rarefaction), a contact discontinuity (at 
X ~ 0.475 when t = 0.55), and three right-going waves (slow shock, Alfven wave, and fast shock). Because of the 
small relative velocity of the different waves emerging from the break-up of the initial discontinuity, which generates 
extremely narrow structures, this test is rather demanding for the different numerical algorithms which resolve the 
wave structure either totally (as it is the case of the FWD solver) or partially (as in the case of HLLC). 

As we show in Fig. [6l all the approximate solvers yield quantitatively similar results across the fastest waves (inde- 
pendent on whether they are heading to the left or to the right), since the numerical flux has to be consistent with 

the upwind update of the state vector, which happens either if A'' > 0, Vp, in which case F — F(Ul), or if < 0, Vp, 

in which case F = ¥{\Jr) (see, e.g.. lAlov et"aIlll999D . 

More obvious differences are noticeable across the contact discontinuity (Fig. [71 left), where both HLLC and FWD 
solvers are closer to the exact solution than the HLL solver. However, because of the extra dissipation produced by 
the rotation of the transverse components of the magnetic field, an unphysical undershoot is produced in the rest-mass 
density (Fig. [3 left). 

The narrow structure generated between the slow shock and the rotational wave (0.7 x < 0.725) is much sharply 
resolved with the HLLC and FWD solvers than with HLL, but the superiority of the Roe-type solution is obvious in 
the middle and right panels of Fig. [71 For the working resolution of 800 uniform zones, none of the three schemes 
represents adequately the extremely narrow structure confined between the left-going Alfven and slow rarefaction 
(0.185 ^ X < 0.19). Indeed, this is to be expected, since re solving this region is a challenge, not only for RMHD 
codes, but also for the exact solver of iGiacomazzo fc Rezzollal (2006) (actually, this is the reason for the impossibility 
of obtaining an exact solution with an accuracy better then 3.4 x 10"'*; c.f.. IGiacomazzo fc Rezzollall2006[ ). 

The overall similar quality of the results obtained for this test using the HLLC and the FWD solvers is also reflected 
in the Ll-norm errors, which are virtually the same in both cases, independent of the numerical resolution employed 
(Fig. [5l upper right). As expected, HLL deviates more from the exact solution and possesses Ll-norm errors ^ 24% 
and ~ 19% larger than those of the FWD and HLLC solvers, respectively. 

* In the following we refer to this combination as the third order scheme 
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For future reference, we provide the ratio of CPU times for this test using the first order scheme, which are tuLL '■ 
tahhc '■ ^FWD = 1 : 1-03 : 3.15. 

8.1.4. Shock Tube 3: Colliding streams 

Figure [5] shows the results of a strong relativistic shock reflection test with two streams approaching each other at 
a ve locity — 0.999. Th e test was proposed by .Balsara (200 1 ,). and it ha s also been used for code validation in, 
e.g. JDel Zanna et all (|2003D ; lLeismann et all (|2005D : iMignone fcBodol ()2006[ ). and MUB09. This setup results in two 
fast and two slow shocks. At the initial collision point (a; = 0.5) a certain amount of "wall heatin g" occurs, which is a 
numerical pathology of approximate Riemann solvers (e.g. INohlll98"7l : IDonat fc Marauinalll996f) . The problem arises 
because an excess of entropy is generated at the collision point at i = 0, which can diffuse numerically only slowly, 
because the fluid is at rest at that point. Higher order reconstruction schemes help to confine the problem to a small 
number of points initially, but generate less diffusion. Hence the "hole" in the rest mass density is deeper with the 
third order scheme than with the first order one. In a close analogy to the HLLD solver of MUB09, the FWD solver 
displays a density undershooting which is smaller than that corresponding to the HLLC solver (the error with respect 
to the analytic solution at a; = 0.5 is 21% and 32% for the FWD and the HLLC solvers, respectively; Fig. [9|). The 
larger numerical diffusion of HLL yields a relative error at the former point of only 10%, but this comes at the cost of 
requiring more zones to resolve both slow shocks. 

The smallest Ll-norm errors correspond systematically to the FWD solver, followed by the HLLC solver and, finally, 
the largest errors are obtained using the HLL solver (Fig. [5] bottom left). In this particular test, the CPU times scale 
as ^HLL : ^HLLC : ipwD = 1 : 1-06 : 2.91. 

8.1.5. Shock Tube 4: Generic Alfven test 

iGiacomazzo fc RezzoUal ()2006D proposed the so-called Generic Alfven test, which MUB09 have also adopted as a 
benchmark for their HLLD solver. The initial discontinuity results into a contact discontinuity which separates a fast 
rarefaction wave (at x ~ 0.05), a rotational wave (at x ~ 0.44) and a slow shock (at x ~ 0.46), from a slow shock (at 
X ~ 0.56), an Alfven wave (at x ~ 0.57) and a fast shock (at x ~ 0.97). 

The exact solution of this test, together with the results for the FWD, HLLC and HLL solvers are shown in Fig. [TOl 
The finest structure in this test, associated with the rotational discontinuities which travel very close to the slow 
shocks, is barely resolved by the first order scheme using the Roe-type solver at the working resolution of 800 uniform 
numerical zones fFig. [TT|) . This is the reason for the considerably smaller Ll-norm errors produced by the FWD solver 
as compared with HLL and HLLC (Fig. [5] bottom right). 

For completeness, we also provide with the CPU times needed to run this test for the different solvers ^hll : ^hllc : 
ipwD — 1 : 1-03 : 3.01, and with the ratio of Ll-norm errors obtained with the HLL, HLLC and FWD solvers, at the 
maximum resolution employed in Fig. [51 which is 1 : 0.66 : 0.38. 

8.1.6. Ultrarelativtstic shell 

With the previous one dimensional tests, we have demonstrated the effectiveness and strength of accurate Riemann 
solvers in the description of complex relativistic magnetized flows, although, as we pointed out in Sect. I8.1TT1 the use 
of high-order schemes tends to reduce the differences between simpler, more diffusive Riemann solvers (like HLL) and 
more elaborate ones (like our FWD solver). However, it remains true that, in multidimensional applications, even in 
combination with high order schemes, it may pay off to use more computationally expensive solvers, which reduce 
substantially the numerical errors without the need of employing huge, perhaps in practice unreachable, numerical 
resolutions. 

In this section we will show that the new Roe-type solver can dodge the difficulties inherent to some realistic 
astrophysical flows, where the other two solvers fail. For this purpose, we have set up a magnetized ultrarelativistic 
shell in spherical coordinates, which moves radially with a Lorentz factor Wq ~ 15. The magnetization parameter is 
(To = &^/47rp = 0.1, corresponding to an azimuthal (toroidal) magnetic field B'^ = 14.56 perpendicular to the purely 
radial velocity (note that in Table 1, where the shell state is denominated M, the vector components {x, y, z) correspond, 
in th is test, to the spheric al components (r, </>, 0)). With this set up we model the ejecta in gamma-ray burst afterglows 
as in iMimica et al.l (|2009f ) , where a pure one-dimensional treatment is adequate because of the ultrarelativistic speed 
of propagation of the ejecta, in spite of the fact that the ejecta is axially symmetric about the direction of propagation, 
but it may be heterogenous in the 0-direction. Therefore, we consider here the evolution of a representative hollow 
wedge of the ejecta located at a small polar angle 6*0 from the spherical axis. The shell has an initial radial thickness 
Aq = 10^*^ cm (Aq = 0.1 in the code units), and moves in a uniform medium whose density is much smaller than that 
of the shell (state R in Table 1). Behind the shell (state L in Table 1), the medium is almost as light as that ahead, 
but moves at the same speed than the overdense shell. 

Because of the numerical difficulty of the proposed set up, we have to run this test with 32000 uniform numerical 
zones in the radial direction, with the third-order scheme and a CFL number 0.2. The results at t = 3.5 (about 
10^ timesteps) can be seen in Fig. [121 just before the HLLC scheme fails to continue the evolution. At this stage 
the breakup of the initial sharp edges of the shell has produced a (reverse-)shock-contact-(forward-)shock structure 
(8.501 X 10^®cm< R < 8.512 x 10^^ cm) at the front edge and a shock-contact-rarefaction at the rear one (8.393 x 
10^^cm< R < 8.413 x 10^® cm)^. We point out that, because of the spherical expansion, the intermediate states 

® Note that, since the magnetic field is perpendicular to the velocity, there are only three waves emerging from the discontinuities -and 
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between the two shocks in the front of the sheU are not uniform. 

The reason for the failure of the HLLC solver is the high frequency noise produced in the rear contact discontinuity, 
which extends with time to the whole rarefied shell (see the gas pressure in the region 8.401 x 10^^ cm < R < 8.413 x 
10^^ cm; Fig. [T^ . This noise is practically nonexistent when we use the FWD solver. Another difference is noticeable 
at the reverse shock (i? ~ 8.501 x 10^^ cm), where the unphysical undershooting is deeper when using the HLLC solver 
than when the FWD solver is employed. 

The test we have presented here argues in favor of the use of more elaborate (although more expensive computa- 
tionally) Riemann solvers for ultrarelativistic applications. 

8.2. Two-dimensional tests 

We employ the framework provided by the RGENESIS code (iLeismann "erall[2f)n5h to set up a pair of two dimensional 
test problems, with the aim of illustrating the multidimensional performance of our new Roe-type solver and the 
proposed renormalised full wave spectral decomposition. For the multidimensional tests we use the third-order scheme 
and the constraint transport method to update the magnetic field flux. Although there is no analytic solution to 
verify against, our results can be compared with those produced by other authors. As these tests encompass all the 
degeneracies of the RMHD eigensystem, they are a challenge for our linearized, Roe-type Riemann solver. 

8.2.1. Cylindrical explosion 

The cylindrical explosion test consists of a strong shock propagating into a magnet ically domina t ed medium. Result s 
for different cylind r ical explosion problems i n RMHD have been p u blished by, e.K.,lDubail (19911 ): Ivan Putt"enl l)1995[ ): 



iKomissarovl (|1999[ ) ; iDel Zanna et al.' (!2003D ; ILeismann et all (|2005D ; iMignone fc Bodol (|2006h . We have chosen a setup 
very similar to that in [Komissarov (1999): a cylinder of high pressure and density is located in the center of a square 
Cartesian grid, which initially contains a uniform, strong magnetic field. Simulations have been carried out on grids 
of 400 X 400 and 800 x 800 zones, spanning a region of 12 x 12 units of area. In the center of the grid there is a circle 
of radius 0.8 where p = 10~^ and p = 1. Between a radius of 0.8 and 1.0 the values exponentially decrease to those 
of the homogeneous ambient medium {p = 10"'' and p = 3 x 10~^ ). Initially, the magnetic field is i^a, = 1, and the 
velocity is zero everywhere. We have c hosen this pa r ticula r explosion set up among several other possible because it 
is the strongest cylindrical explosion in Komissarovl (|1999() . Several other possible initial states have computed with 
the new Roe- type solver in Anton (2008). 

The test runs with different solvers have been performed using our third-order scheme and a Courant numb er 0.4, 
for th e highest resolution runs, and 0.6 f or the lower resolu tion ones. To compute the most extreme of the Komissarovl 
(|1999D tests, the energy fix proposed by IMignone fc Bodol (2006) (cf., their equation 76) has been used. 

The explosion of the central overpressured region is notably confined by the presence of the strong magnetic field, 
which hinders the expansion of the fluid in the y— direction, and yields a pair of twin jets propagating along the 
x— direction. The outer fast shock has an almost spherical shape, with a hump along the x-direction, where the 
maximum Lorentz factor of the expanding fluid is reached^^. This maximum Lorentz factor attained by the models at 
t = 4: depends on the numerical resolution and on the solver, and ranges from W,nax — 3.5 to 4.03, with HLL yielding 
the smaller values of Wmax and the FWD solver the largest ones. On the other hand, at higher resolution Wmax is 
about a 10% larger than at low resolutions. Figure [Ql shows that the FWD solver produces the sharpest and best 
resolved features in the distribution of the rest-mass density at t = 4. However, the results we obtain using HLLC 
and FWD solvers are very similar, particularly at high resolution, the only noticeable differences being in the narrow 
horizontal region of high density (along which the explosion happens) and the shape of the bow shocks. Compared 
with HLL, the FWD solver delivers sharper profiles in the discontinuities even with half the resolution (compare in 
Fig. [Impanels c and d). 

In this test, whose dynamics is dominated by the fast shocks emerging from the break-up of the initial discontinuity 
the ratios of CPU times are tuLL '■ ^hllc ■ tpwD = 1 : 1-07 : 1.52 and ^hll : ^hllc : ^fwd = 1 : 1-07 : 1.86 at resolutions 
of 400 X 400 and 800 x 800, respectively. This means that the FWD solver is less demanding computationally than in 
all the ID cases analyzed in Sect. 18.11 The advantage of using the FWD solver over the HLLC one is questionable, 
provided the similitude of their respective results. However, it is possible to find a trade-off between the numerical 
efficiency that the HLL solver yields and the accuracy that the Roe-type solver lends. 

8.2.2. The rotor problem 

IDel Zanna et all (j2003l ) adapted the rotor problem of classical MHD shown by iBalsara fc SpiceJ ()1999l ) and iTothI 
(pooofT to the relativistic case. In a domain [—0.5, 0.5] x [—0.5, 0.5] in the a:y-plane, a circular region of radius Tc = 0.1 
rotates rigidly with an angular velocity lOc = 9.95, i.e., with linear velocities (v^,v^) = ijjc{—y,x). The central circle 
is homogeneous, with a density pc = 10. The medium surrounding the central circle is also homogeneous and has a 
density p = \. The thermal pressure, the magnetic field and the adiabatic index are constant everywhere, pc = 1, 
(i?^i?^) = (1,0), 7 = 5/3. 

We carry our computations until t — 0.4 at resolutions of 128^, 256^, 512^ and 1024^ employing HLL, HLLC and 
FWD solvers. We use the third-order algorithm with a Courant number 0.6. We point out that this test has a 

not seven as in non-degenerate RMHD-. 

One could associate these humps with the bow shocks of both relativistic, oppositely-moving, twin jets. 
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resolution-dependent complexity since the maximum Lorentz factor in the set up of the initial model depends on the 
resolution. Theoretically, the maximum Lorentz factor in the initial set up is reached at the circumference that bounds 
the central circular region (i.e., at r = rc), where W^m^x = (1 ~ r^w^)"^/^. The finer the resolution, the closer is the 
maximum Lorentz factor in the grid to since the center of the cells crossed by the circle r = rc is closer to Vc- 

Even at the coarsest resolution employed (128^), both solvers, HLLC (Fig. [Tib ) and FWD (Fig. [TiH). capture 
adequately the two rotational discontinuities that are generated from the points (x, y) = (0, ±0.1), where the magnetic 
field is initially parallel to the velocity field. The HLL solver smooths out such discontinuities, even at the highest 
resolution we have considered (1024^; Fig. [T4k). 

Both FWD and HLLC solvers yield very similar results independent of the numerical resolution, but HLLC is 
between a 40% and a 70% faster at high resolution. We note that increasing the resolution brings a reduced ratio 
^fwd/^hll (Table 2), which implies that it is more favorable to use a more accurate solver (FWD) at moderate or 
high resolution than a more diffusive one (HLL). 

9. SUMMARY AND CONCLUSIONS 

In the first part of the paper, we present renormalized sets of right and left eigenvectors of the flux vector Jacobians 
of the RMHD equations, which are regular and span a complete basis in any physical state including degenerate ones. 
Our starting point are the expressions of right and left eigenvectors in covariant variables (Sects. 13.31 and 13. 4p . The 
renormalization procedure relies in the characterization of the degeneracy types in terms of the normal and tangential 
components of the magnetic field in the fluid rest frame (Sect. I3.2p . Type I (II) degeneracy occurs whenever the 
magnetic field is purely tangential (normal) to the Alfven wave front. Then, proper expressions of the renormalized 
right and left eigenvectors in conserved variables are obtained through the corresponding matrix transformations 
(Sects. [5] and [6|). Our work completes previous analysis (Balsara 2001, Komissarov 1999, Koldoba et al. 2002) that 
present different sets of right eigenvectors for non-degenerate and degenerate states, and can be regarded as a relativistic 
generalization of the work performed by Brio & Wu (1988) in classical MHD. Our theoretical analysis is completed 
with a brief note about the non-convex character of the relativistic MHD (Sect. [4]). 

The second part of the paper deals with the development of a linearized Riemann solver based on the renormalized 
spectral decomposition of the equations just obtained. Here our work departs from previous works along the same line 
(Balsara 2001, Komissarov 1999, Koldoba et al. 2002) in that we build the numerical fluxes from the (renormalized) 
analytic right and left eigenvectors in conserved variables. Technical details of the numerical code are given in Sect. [71 
Intensive testing against one- and two-dimensional numerical problems (see Sect. [5]) allows us to conclude that our 
solver is very robust. When compared with a family of simpler solvers (HLL, HLLC, HLLD) that avoid the knowledge 
of the full characteristic structure of the equations in the computation of the numerical fluxes, our solver turned out 
to be less diffusive than HLL and HLLC. This is clearly seen in the tests presented in Sect. 18.1.11 with both first 
and third order algorithms, that involve the description of isolated contact and Alfven discontinuities. This is an 
expected result since Roe- type Riemann solvers are based in the full wave decomposition of Riemann problems. This 
conclusion is confirmed by the numerical results shown through Sects. l8.1.2j[S7l.5l Comparison with HLLD solver is 
not as favorable. Here, the expected superiority of Roe-type Riemann solvers with respect to central TVD schemes 
(see, e.g., the discussion in Toro 1997) is put into question since the increased accuracy of the five- wave HLLD solver 
on one hand, and the large amount of operations involved in the computation of numerical fluxes in the FWD solver on 
the other hand, lead to comparable numerical accuracies. This large amount of operations needed by the FWD solver 
makes it, in principle, less efficient computationally than those of the HLL family. However the relative efficiency of the 
FWD solver increases in multidimensional simulations since the number of zones needed to achieve a given accuracy 
is much smaller in the case of the FWD solver than in that of more diffusive solvers. 

It is well know, that choosing the adequate restrictions on characteristic variables serves optimally the purpose of 
setting boundary conditions. The computation of characteristic variables needs of the knowledge of the complete set 
of left and right eigenvectors. Thus, in spite of the efficiency and complexity arguments brandished in the previous 
paragraph, we point out that the set of computed left and right eigenvectors has the potentiallity of being used 
to properly set boundary conditions in a number of numerical applications. Work along this line will be presented 
elsewhere. 
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APPENDIX 



CHARACTERISTIC WAVESPEEDS DIAGRAMS AND DEGENERACIES 



Contrarily to what happens in the case of (classical and relativistic) hydrodynamics, in a homogeneous magnetofluid 
at rest the presence of a magnetic field leads to a dependence on direction in the propagation speed of waves. This 
dependence can be visualized by mean s of the characteristic wavespeeds diagram, equivalent to the phase speed 
diagrams introduced by Friedrichs (see peffrev &: Taniutlll964[ )). These diagrams show the normal speed of planar 
wavefronts propagating in different directions, the speed is given by the distance between the origin and the normal 
speed surface along the corresponding direction. 

In a recent paper, Keppens and Meliani (2008) have revisited the theory for linear RMHD wave propagation showing 
the equivalence with the characteristic speed approach, and the connection between the phase and group speed diagrams 
by means of a Huygens construction, in arbitrary reference frames. In this Appendix, we discuss our results obtained 
in the representation of the characteristic speeds for fluids under different thermodynamical conditions, focussing in 
the layout of degeneracies as a function of the fluid's motion. Our results can be directly compared with those shown 
in Keppens and Meliani (2008). 

Top-left panel of Fig. [15] displays a cut in the a;?;-plane of the characteristic wavespeed surfaces corresponding to a 
homogeneous magnetized fluid at rest. In this case, the surfaces exhibit two symmetries. The first one is a rotation 
symmetry around the direction of the magnetic field, which justifies our choice of displaying a plane containing the 
magnetic field (along the x-axis in the figure). The second one is a mirror symmetry across the orthogonal plane to 
the magnetic field. The two types of degeneracy present in classical and relativistic MHD are found on this orthogonal 
plane (Type I) and along the magnetic field direction (Type II). Top-left panel of Fig. [15] shows a particular subcase 
of Type II degeneracy, namely the one in which the speeds of the Alfven waves coincide with the corresponding fast 
magnetosonic waves. Top-right and bottom-left panels of the same figure display the remaining two subcases (see 
the corresponding figure captions for the fluid conditions). In the top-right panel the speeds of Alfven waves along 
the magnetic field direction coincide with those of slow magnetosonic waves. The triple degeneration case in which 
Alfven, slow and fast magnetosonic waves propagate at the same speed along the magnetic field direction is shown in 
the bottom-left panel. 

For fluids in motion, the diagrams loose the symmetry properties described in the previous paragraph. In the case 
of motion along the direction of the magnetic field, as in the bottom-right panel of Fig. [13 the diagram still displays 
the rotation symmetry around the magnetic field direction and the degeneracies of types I and II still occur in the 
orthogonal plane to the magnetic field and along the magnetic field, respectively. When the fiuid velocity is not aligned 
with the magnetic field, the rotation symmetry around the magnetic field direction is lost. Figure [16] displays cuts of 
the characteristic wavespeeds surfaces on the xy and xz planes for a magnetized fluid under suitable thermodynamical 
conditions moving along the y-axis (the magnetic field is still along the x axis). In this particular case in which the 
motion is orthogonal to the magnetic field direction, the mirror symmetry across the yz-plane is retained. Type I 
degeneracy is also found on this plane. Concerning Type II degeneracy, it is remarkable to note that the degeneracy is 
not present in the a;z-plane, whereas in the xy-plane appears along two different directions. Moreover, contrarily to the 
classical MHD case, along each of these directions only one of the Alfven waves is degenerate with the corresponding 
(fast) magnetosonic wave. 

In the most general case (Fig. I17|) a mirror symmetry across the plane defined by the magnetic field and fiuid velocity 
vectors (the xy-plane in these figures) is still retained. Again, Type II degeneracy appears along two different directions 
and only one of the Alfven waves is affected by the degeneracy. 

PROOF OF THE PROPORTIONALITY BETWEEN THE MAGNETOSONIC RIGHT EIGENVECTOR AND THE MODULUS 

OF THE TANGENTIAL MAGNETIC FIELD 

In this appendix we will proof that the components of the magnetosonic eigenvector, equation (|69|) . are proportional 
to \bt\. For a magnetosonic wave, and similarly to the Alfven wave case, we can decompose the magnetic field in 
perpendicular and parallel part to the magnetosonic front. If we do so, 5" is written as 



and, from here. 



Substituting this decomposition in the quartic equation. Mi = 0, it is straightforward to obtain the following relation. 



6" = 



B 



(0" + au") + b\ 



(Bl) 



G + a2 




(B3) 



Let us now transform the components of the magnetosonic eigenvectors (equation [69]) using the above expressions. 
For the first four components we have 
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Similarly, for the following four components 



And for the ninth component of the eigenvector, 



l + (B5) 



Since the tangential magnetic field can be written as the product of its modulus and a unitary vector, the 
proportionality is proven. 
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Fig. 1. — Density and y component of magnetic field from an isolated contact (left panel) and a rotational (right panel) wave at t = 1 (as 
in MUB09). The different symbols show results computed with the FWD solver (triangles), the HLLC solver (diamonds) and the simpler 
HLL solver (plus signs). Note that both FWD and HLLC solvers are able to capture exactly an isolated contact wave, keeping it perfectly 
sharp without producing any grid diffusion effect. HLLC can capture the contact wave but not the rotational discontinuity, whereas HLL 
spreads both of them on several grid zones. The Roe-type solver captures with very small diffusion the rotational wave. 
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Fig. 2. — Same as in Fig. [T]but using a third order accurate Runge-Kutta time integrator and PPM spatial reconstruction. The diffusion 
of the HLL solver across the contact wave is drastically reduced. The same comment applies for both HLLC and HLL across the rotational 
wave. We note the small oscillations behind the contact (in the HLL and HLLC cases) resulting from the higher order of the method. 
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Fig. 3. — Rclativistic Brio-Wu shock tube test at t = 0.4. Computations are carried on 400 zones using the FWD (dashed-dotted Unc), 
HLLC (dashed hne), HLL (dotted-crossed Une) and exact (sohd hne) Riemann solver, respectively. The top panel shows, from left to 
right, the rest-mass density, gas pressure, total pressure. In the bottom panel the x and y components of velocity and the y component 
of magnetic field are shown. The compound wave in the numerical solution between the fast left-propagating rarefaction and the contact 
discontinuity (at x « 0.51 in the plots) is not captured in the analytical solution provided by Giacomaazo & RezzoUa (2006) that avoids 
this kind of solution by construction. 
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Fig. 4. — Zoom of the central region of Fig[3l focusing on the compound wave (at x ~ 0.51), the contact discontinuity (at x ~ 0.60) and 
the right-going slow shock (at x ~ 0.65). From left to right, we show the rest-mass density and the two components of velocity. 




Fig. 5. — L-1 norm errors (in 10^) for the four shock tube problems of Table 1 as function of the grid resolution A.x = 1/Nx. The different 
combinations of lines and symbols refer to FWD (solid, diamonds) , HLLC (dashed, crosses) and HLL (dotted, plus signs) . 
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Fig. 6. — Second shock tube problem. From left to right, the top panel shows rest-mass density, gas and total pressure. The middle 
panel displays the three components of velocity. The bottom panel shows the Lorentz factor and the transverse components of magnetic 
field. Solid, dashed-dotted, dashed and dotted lines are used to identify results computed with the exact, FWD, HLLC and HLL solvers, 
respectively. 
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Fig. 7. — Left panel: closer look of the rest-mass density of Fig. |6] around the contact wave. Middle and right panels; enlargements of 
the z component of velocity and y component of magnetic field around the right-going slow shock and Alfven discontinuity, respectively. 
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Fig. 8. — Results of the shock tube 3 (Table 1) using different solvers. From top to bottom, left to right, the panels show rest-mass 
density, gas pressure, total pressure, x and y components of velocity and y component of magnetic field. The exact solution of this Riemann 
problem is plotted with a solid line. Dash-dotted, dashed and dotted lines refer to the results obtained with FWD, HLLC and HLL solvers, 
respectively. 
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Fig. 9. — Closer view of the central region of Fig|H] Type lines have the same meaning as in Fig[7] The rest-mass density is represented 
in the left panel, where the wall heating pathology is notorious. Central and right panels show the profiles transverse components of the 
velocity and of the magnetic field, respectively. 
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Fig. 10. — Same as Fig. |6] but for the general Alfven test (Table 1). 
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Fig. 11. — Enlargement of the central region of Fig llOl where the profiles of the rest-mass density (left), (middle) and are displayed. 
The contact wave at a; ~ 0.52 separating the two slow shocks is better captured by the FWD solver. Both HLLC and FWD solvers 
undershoot the true solution to the right of the contact wave. The Alfven modes at x ^ 0.44 and x 0.58 are not fully smeared at this 
resolution using the FWD solver, although they are obviously unresolved using 800 uniform numerical zones. 
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Fig. 12. — Ultrarelativistic shell at 4 = 3.5 (simulation units, which correspond to 1.1 X 106s) run with 32000 numerical zones (Table 1). 
The left and right panels show the results using HLLC and FWD solvers, respectively. Different line styles are used to display a number 
of physical variables. 
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Fig. 13. — The rest-mass density distribution of the two-dimensional cyhndrical explosion problem at t = 4.0. The upper (bottom) panels 
have been run with a grid of 800 X 800 (400 X 400) numerical zones, a CFL number 0.4 (0.6), PPM spatial reconstruction and a third-order 
accurate Runge-Kutta integration scheme. Left, middle and right panels show the results obtained using FWD, HLLC and HLL solvers, 
respectively. 




Fig. 14. — The rest-mass density distribution of the two-dimensional rotor problem at t = 4.0. The upper (bottom) panels have been run 
with a resolution of 1024^ (128^). Left, middle and right panels show the results obtained using FWD, HLLC and HLL solvers, respectively. 
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Fig. 15. — Characteristic wavespeeds. Ideal gas with 7 = 4/3. Dashed lines correpond to fast magnetosonic waves, continuous lines to 
Alfven wavespeeds, and dash-dotted lines to slow magnetosonic waves. Top-left panel: Fluid at rest with p = 1.0, e = 1.0, B^' = 5.0, 
By = 0.0, = 0.0. Top-right panel: Fluid at rest with p = 1.0, e = 50.0, = 5.0, Bf = 0.0, = 0.0. Bottom-left panel: Fluid at 
rest with p = 1.0, s = 37.864, B^ = 5.0, B" = 0.0, B^ = 0.0. In these three panels, the curve associated to the entropy wave degenerates 
in a point located at the origin, and Type I degeneracy is along y axis, whereas the three subcases of Type II degeneracy are along the 
X axis. Bottom-right panel: Fluid state as in the top-right panel but moving along the x-axis at a speed = 0.9. Entropic wavespeed 
corresponds to dash-triple-dotted line. Type I and II degeneracies are again along y and x axis, respectively. 
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Fig. 16.— Characteristic wavespeeds. Ideal gas with 7 = 4/3, p = 1.0, e = 1.0, i)^ = 0.0, = -0.50, -u^ = 0.0, = 5.0, = 0.0, 
= 0.0. Line types as in Fig. 1151 Left panel: xy-plane. Type I degeneracy is along y axis. Type II degeneracy is along the directions 
pointed by the arrows. Right panel: xz-plane. Type I degeneracy is along z axis. There is no Type II degeneracy on this plane. 
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Fig. 17.— Characteristic wavespeeds. Ideal gas with 7 = 4/3, p = 1.0, e = 1.0, v"^ = 0.634, = 0.634, = 0.0, = 5.0, BV = 0.0, 
B^ = 0.0. Line types as in Fig. 1151 Left panel: xy-plane. Type I degeneracy is along y axis. Type II degeneracy is along the directions 
pointed by the arrows. Right panel: xz-plane. Type I degeneracy is along z axis. There is no Type II degeneracy on this plane. 
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TABLE 1 

Initial conditions for, the test problems of Sect. 18.11 



Test 


State 


P 


Pg 




vy 






BV 




Time 


Zones 


Contact Wave 


L 
R 


10 
1 


1 
1 






0.7 

0.7 


0.2 
0.2 


5 
5 


1 
1 


0.5 

0.5 


1 


40 


Rotational Wave 


L 

r> 
irv. 


1 

1 
± 


1 

T 
1 


0.4 

U.O / / O'^ 1 


-0.3 
— u.^ozooy 


0.5 

n AOAA on 


2.4 

9 A 


1 

n 1 

— U. i 


-1.6 


1 


40 


Shock Tube 1 


L 

R 


1 

0.125 


1 

0.1 














0.5 
0.5 


1 

-1 






0.4 


400 


Shock Tube 2 


L 
R 


1.08 
1 


0.95 

1 


0.4 
-0.45 


0.3 
-0.2 


0.2 
0.2 


2 
2 


0.3 
-0.7 


0.3 
0.5 


0.55 


800 


Shock Tube 3 


L 

R 


1 
1 


0.1 
0.1 


0.999 
-0.999 










10 
10 


7 

-7 


7 

-7 


0.4 


400 


Shock Tube 4 


L 
R 


1 

0.9 


5 

5.3 






0.3 




0.4 



1 
1 


6 
5 


2 
2 


0.5 


800 


UR shell 


L 
M 
R 


lO-'^ 
1 

3.297 X 10-* 


10-'' 
2.5 X 10-^ 
8.2424 X lO-'^ 




0.997775 
0.997775 



















4.60706 








3.5 


32000 



Note. — The last two columns give, respectively, the time at which the solution is shown in the figures and the number of numerical 
zones used in the computation. The last test is set up in spherical coordinates, and the vector components {x. y, z) correspond to the 
spherical ones (r, (p, 9) 
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TABLE 2 

CPU TIMES FOR THE TWO-DIMENSIONAL 
ROTOR PROBLEM. 



Resolution tnLLc/tHLL <fwdAhll 

128^ 1.75 2.53 

256^ 1.45 2.84 

512^ 1.08 1.50 

1024^ 1.13 1.93 



Note. — Execution times for test cases 
run with the HLLC and FWD solvers are 
normalized to those corresponding to the 
HLL solver. In all cases the third-order 
scheme with a Courant number 0.3 has been 
used 



